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18. SUPPLEMENTARY NOTES (Continued) 


EDITORS’ NOTE 


Volumes 13 to 17 were originally published by SAI to describe the atmo- 
spheric, geomagnetic, and high-altitude energy deposition and neutral 
heave models for ROSCOE. This whole section of code, when associated 
with an appropriate DRIVER subroutine, operated as a package that ran 
independently of the rest of the ROSCOE structure. Provision was also 
made, within this high-altitude package, for two completely independent 
descriptions of atmospheric heave, each with its own description of 
atmospheric chemistry. 


When GRC incorporated this section of code within the ROSCOE frame- 
work, some modifications were necessary, which means that some of 

the descriptions in Volumes 13 to 17 are inappropriate to ROSCOE as 

it now exists. In particular, the NRL heave routines (deck NRLHYD) 

and associated chemistry (deck NRLCHM) are not presently used in 
ROSCOE. Three other subroutines are different: subroutines ATMOSU, 
EIF, and XTCOEF correspond to the ROSCOE subroutines ATMOS, 
EXPINT, and WDXP, respectively. With these exceptions, the subroutines 
described in Volumes 13 to 17 correspond exactly to those currently in 
ROSCOE. 


20. ABSTRACT (Continued) 


method of either difference equations or differential quadrature, 
respectively. Both models have automatic rezone capability. The 
geocentric columns defining the geometry for the calculations are : 
described in terms of an arbitrarily positioned and oriented quadrupole 
coordinate system. Each Lagrangian cell or point is characterized by 
not only the hydrodynamic properties but also a set of chemistry 
quantities. The chemistry is loosely coupled to the hydrodynamics. 
Herein are presented details of the quadrupole coordinate system, the 
working form of the hydrodynamic equations and their initialization and 
methods for solution, the results of a test problem for a large-yield 
event at 200-km altitude producing motion in a linear array of six 
columns. 
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1. INTRODUCTION 


tod DEFINITIONS 


For the purpose of the New Code development, the phrase at- 
mospheric heave, or simply heave, refers to the motion (the general 
upwelling) of the atmosphere above approximately 100-km altitude, 
induced by energy deposition. Asa first approximation, the modeling 
of the overall phenomenon has been broken into two portions: the neutral 
component (Model No. 15a) and the ionized component (Model No. 15b). 
These two portions are, of course, closely related, but they have been 
partially decoupled in this early modeling effort. It was the Commit- 
tee's plan |see Vol. 13] to model the neutral component (SAI's re- 
sponsibility) without attention to the ionized component (PDI's responsi- 
bility). The modeling of the ionized component involves a structure that 
uses the results from the neutral-component motion. It should be noted 
that the ionized component being considered is a weakly-ionized com- 
ponent in the UV halo region, to be used as a basis for describing stri- 
ations (Model No. 17), and not that in the fireball region which will be 


described in still another model, No. 1la. 


1.2 HISTORICAL BACKGROUND (BRIE F) 


La 2s2 Early Developments in Heave 


Atmospheric heave is, perhaps, the single most important 
multiburst effect. The phenomenon was discovered in 1966, first de- 
scribed by R. W. Lowen in a meeting held at the Rand Corporation in 


May 1966 and reported in HL-66a. (The phrase, atmospheric heave, 


came into general usage after Hendrick introduced it |He-67a].) The 


«+ Ore 


first calculations were done using a hydrostatic model {HL-66a]|, but it 
was recognized early that this neglected a number of effects that were 
likely to be important. Consequently, better calculations were next 
performed using a one-dimensional Eulerian hydrodynamics code 
|HL-67b, HL-67d, HL-69b]. These calculations gave an idea of the 
time scale for the effect to develop and showed two additional features: 
a Strong upward-moving shock developed early, and an overshoot oc- 


curred in the density profile. 


Because it was recognized that two-dimensional (or, ina mul- 
tiburst environment, three-dimensional) effects might well matter, 
more elaborate calculations in two dimensions were performed with the 
SHELL code in collaboration with personnel at AFWL, firs* for x rays 
alone {Hamlin et al., First Meeting of RANC Working Group held at 
G. E. TEMPO, Santa Barbara, 1-4 October 1968; HL-69b, Section 2: 
Comparison of One- and Two-Dimensional Hydrodynamics Calculations 
of X-Ray~Produced Atmospheric Heave; BF-69], and later with x-rays, 
UV, and debris included |HL-69c, -69d, -69e|. Several important con- 
clusions were apparent: soft UV was important |Ha-69b]|, radial motion 
was important and was the principal cause of later subsidence, and the 
shock and overshoot shown in the one-dimensional calculations were not 


present. 


Bs bare RANC IV Heave Models 


For systems codes, heave was first incorporated into RANC 
IV (the formal development of which began in October 1968 with the first 
of about a half-dozen working-group meetings that occurred every sev- 
eral months). The basic portion of the heave model ultimately devel- 


oped for RANC IV was a set of equations which described the trace- 


particle trajectories obtained from a series of one~dimensional, 


vertical-column hydrodynamic-code calculations. Multiburst effects 
were estimated by summing the air-parcel displacements caused by 
each of the bursts; the displacement increment given by the nth burst 

to a Specifie air parcel depended on the altitude of that parcel at the 
time of the nth burst. The heave model also had a procedure to approxi- 
mate the multiburst heave-reduction effect, in which the displacement 


caused by a later burst is reduced because of the extra mass lofted high 


in the atmosphere by earlier bursts |MF-73a 


. 


Two.serious deficiencies in the original heave model were 
recognized |MF-73a|: slow computational speed and discontinuity of 
motion across fireball surfaces. The number of multiplications required 
for a multiburst heave calculation varied as the square of the number of 
bursts for a small number of bursts and as the cube for a large number. 
Lack of close coupling between the heave and fireball models was mani- 
fested by the ghost-fireball phenomenon, which arose because the heave 
model did not predict as high a rise rate for air near the fireball as did 


the fireball motion model |MF-73a |. 


Fischer |Fi-71c] considered some possibilities for an im- 


proved heave model for RANC IV. 


A vastly more efficient heave model for RANC IV has also been 


described by McNamara and Finn |MF-73a |. 


1 a. Old-ROSCOE Heave Model 


For the old- ROSCOE code, it was the intent to find a structure 
for modeling heave which would (1) reproduce the essential predictions 
of the large three-dimensional two-fluid MHD code calculations (which 
would form the data base for heave) and (2) run reasonably fast ona 
computer |YS-73, ZS-73a|. The general approach was to use a quasi- 


hydrostatic model with faked-in time dependence through adjustable 


eR Re 


ae RE ae 


parameters, values of which were selected so that the density profiles 
approximated those found from two- and three-dimensional hydrody~ 
namics codes. Results from a preliminary version of the model, with 
comparisons to hydrodynamics-code results, were presented in April 
1973 |ZS-73a|. The model was in almost constant revision until work 

on the old- ROSCOE code was terminated in November 1973. For exam- 
ple, an improved version of the model, which requires about 1/5 of the 
number of storage quantities as the earlier version |ZS-73a], has been 
described by Stuart [St-73a]. For the latest version see the old- ROSCOE 
final report |LL-75]. 


1.2.4 HARC Heave Models 


For the HARC code, there was less documentation available 
than for the old-ROSCOE code. The statements that follow are based 
on notes taken by D. A. Hamlin on 29 January 1974 during a visit by 
the New Code committee to MRC where D. H. Sowle described the 
HARC heave model. 


In HARC, the heave model describes mass, species, and stri- 
ations (i.e., the gross-scale variety, not those due to fireballs or beta 
tubes). There are three ways to get heave results. The first two ways 
involve interpolations in space and time; the third way is for one event, 
outside the normal battle mesh, for which the treatment is RANC-like 
in nature in that it "recreates the world."’ The third method is always 
in the code; either the first or the second or both are in, at the user's 


option. 


The first method, developed for BTL, uses pre-computed 
phenomenology stored on tapes. Thus, the user must accept what is on 


the tapes, aithough he can override the fireball and can add debris, 


12 
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beta rays, etc. An interface routine can fake some quantities if they 
are not on the tape. Some advantages of pre-computed phenomenology 
are that: 

1. There is none of the delay that normally occurs between 


the time of a new description of phenomenology and the 
time that the results are modeled. 


2. It is the most accurate, in that it obviates the loss of 
fidelity in modeling; thus it can be used to verify models 
used for other purposes. 


3. It is the fastest method, provided the tapes are indeed 
available. 
Some disadvantages of pre-computed phenomenology, which precluded 
it being the only approach to describing heave, are that: 
1. It is non-interactive (i.e. , choices of burst conditions 


must be made ahead of the battle and cannot be made during 
the battle). 


2. It limits the systems analyst to what is on the tapes. 


It requires large storage. 


The second method is described as on-line grid, in which 
multi-one-dimensional Lagrangian hydrodynamics is used. It employs 
a spherical earth-centered coordinate system. The grid contains 
10 x 10 * 13 points, corresponding to approximately 100 km x 100 km 
horizontal resolution and Ah vertical resolution, with Ah varying from 
7 km at the bottom, 93-km altitude, to 100 km at 400-km altitude. The 
mesh resolution is limited by storage and not by running time. A heave 
calculation consists of using the (explicit) Lagrangian hydrodynamics 
plus Scheibe's chemistry modeling (which feeds back energy to the 
hydrodynamics) plus Stagat's striation modeling [described ina BTL 
report and in SC-73a; also see KS-73|. The Lagrangian hydrodynamics 


is taken [partly| from Sappenfield's version of RADFLO in spherical 


geometry; the Richtmyer-von Neumann artificial viscosity is used in 
describing shocks. Time steps are governed by the White instead of 
the Courant condition. The bottom cell is fixed; the top cel) moves 
against a constant pressure, so that the force increases as the area in- 
creases. There is no rezone capability, which would be desirable. 

Two advantages of the on-line grid are that (1) there is no modeling re- 
quired and (2) it is fast. For example, fora 300-sec real-time 
problem, only 10 sec of CPU time on the CDC 7600 (or 50 sec on CDC 
6600 or 130 sec on CDC 6400) is required; the breakdown is about 3 sec 
for hvdrodynamics, 3 sec for chemistry, 2 sec for striations, and 

1 sec for plotting. Results are stored at two times ty and ty ty + A, 
where O=1,..., 20, 30 sec, depending on the time since the last 
burst; information at intermediate times is obtained by interpolation. 


A flag is switchea instead of data being switched in storage when data 


at a new ty are obtained. 


quasi-two-dimensional; it is a fit to two-dimensional hydrodynamics. 
Scheibe's chemistry and Stagat's striation models are again used. 
Elimination of this third method would be desirable but is prevented by 


lack of storage for a larger on-line grid. 


1 ap ea NRL Heave Model 


Subsequent to the initial selection of a heave model as reported 
in Section 1.4 below, D. A. Hamlin and W. S. Knapp visited NRL on 
5 April 1974 where T. P. Coffey described a new type of heave model 
being developed by J. Gardner of NRL. This model, based on the 
differential-quadrature technique described by Bellman et al. | BK-72b], 
had been implemented in a one-dimensional version and was being 


worked upon in a two-dimensional version. 


. 

On 22 April 1974, the COR for the New-Code development, 

P. B. Fleming, requested SAI to evaluate this promising new model, 
noting the following appealing aspects: 

a. The ease of possible rezoning to keep definition in the 
battle space. 

b. The possibility of advancing each vertical array where 
desired at a radar rate, thus eliminating the need for 
storing two arrays at different times. 

c. The fact that placement of these vertical arrays can be 
arbitrarily set up with (perhaps) no need to maintain 
definite separation distances. 

As a result, over the next six-month period, T. P. Coffey 
kindly made available to SAI a series of evolving versions of the NRL 
one-dimensional heave model, including both decks (NRLHV1, 13 May 
1974; NRLHV2, 18 June 1974; NRLHV3, 10 September 1974; NRLHV4, 
13 November 1974) and preliminary memoranda |Ga-74a, Ga-74b, 
GP-74] and a final report {[GP-75]. All versions except the first in- 
cluded the NRL 12-species E- and F-region chemistry package (the 
NRL SIMPLE code |Al-73]). During the period from May to October 
1974, we exercised these codes and made timing runs, reported in the 
bimonthly progress reports to GRC. 

| On 16 October 1974, P. B. Fleming requested SAI to incor- 
porate the NRL heave code into the SAI high-altitude grid module as an 
alternative heave model. We undertook this additional task as a parallel 
effort while we completed the high-altitude grid module containing the 
more conventional difference-equation treatment of the Lagrangian hy- 
drodynamics. Also in October 1974, NRL was requested to incorporate 


an automatic rezone capability into their code, which was included in 


the version NRLHV4 delivered on 13 November 1974. 


jp eas Other Heave Models 


Treatments of heave in WOE II and OPTIR are described in 
SD-73b (Section 5) and MS-73c (Chapter V). 


1.3 RELATION OF HEAVE MODEL TO OTHER ROSCOE MODELS 


The heave mode] will play a central role inthe code. Tenta- 
tively, the heave model (No. 15a) will be affected by models for the 
ambient atmosphere (No. 1), radiation transport and deposition (No. 4), 
possibly low-altitude fireball motion (No. 5a) for fireballs rising into 
the heave mesh, high altitude fireball motion (No. 11a), and chemistry 
(No. 11b and 16). The models that are (tentatively) affected by the 
heave model are those for radiation transport and deposition (No. 4), 
high-altitude fireball properties (No. 11a), high-altitude debris motion 
(No. 12a and 12b), UV generation and deposition (No. 14), heave (ion 
motion, No. 15b), high-altitude chemistry (No. 1b), and striations 
(No. 17). 


1.4 INITIAL SELECTION OF A HEAVE MODEL 


In early 1974 it seemed that the consensus of the New Code 
committee was that each existing heave model was far from ideal, but 
that a model patterned after the HARC on-line grid model would have 
the best chance of becoming operational within the allowed time and re- 
sources. Accordingly, SAI initially put most of its heave~modeling 
effort into a thorough review, revision, and extension of the HARC on- 
line grid model for heave. This choice largely dictated the overall 


structure of the high-altitude grid module. 


2. GEOMETRY OF THE HIGH-ALTITUDE GRID (HAG) MODULE 


In computing the hydrodynamics in a one-dimensional (geo- 
centric) spherical coordinate system, it is not necessary to specify the 
lateral (or angular) extent of the 'column.' However, in depositing the 
energy for an event, which is going to drive the hydrodynamics, it is 
necessary to deal with specifically located boundaries and absolute 
amounts of areal densities and mass. Thus, it is necessary to describe 


the coordinate systems used in the SAI HAG module. 


2.1 QUADRUPOLE COORDINATE SYSTEM 


First, however, we note that the HARC code used an array of 
steradianal columns bounded by meridian planes in the geographic north- 
south direction and by conical surfaces (which intersect the Earth's 
surface in parallels of latitude). The cells ina steradianal column were 
bounded by (geocentric) spherical surfaces. Thus, a HARC cell had 
spherical caps for top and bottom, plane surfaces for east and west 
sides, and conical surfaces for north and south sides. Determination 
of the line segments of a ray path through such cells is extremely 


complex. 


We made a number of improvements and corrections on the 
HARC grid. Two major improvements are (1) replacement of the coni- 
cal surfaces by another set of meridian planes, thereby using what we 
have termed a quadrupole coordinate system (since there are east and 
west poles as well as north and south poles) and (2) replacement of the 
fixed geographic north-south orientation by an arbitrary orientation, so 
that we speak of artificial (or quadrupole) north, south, east, and west 


poles. The quadrupole coordinate system is illustrated in Fig. 1. 
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Geocentric Column in Quadrupole Coordinate System. QN and 
QW are the quadrupole north and quadrupole west poles, and 
O is the quadrupole origin. The quadrupole south and quadru- 
pole east poles are not shown for simplicity. The shaded 
region is a Spherical-surface cross-section through a geo- 
centric column. The quadrupole colatitude (QCZ) and 
quadrupole east longitude (QLZ) of the Point Z are shown. 
The origin need not be at QCZ = QLZ = 90. 


2.2 GEOGRA PHIC AZIMUTH OF QUADRUPOLE SYSTEM 


In practice, we have oriented the quadrupole coordinate system 
so that an artificial north-south meridian will pass through the north 
pole of the magnetic dipole determined by the ambient geomagnetic field 
model |Vol. 3, LS-75]|. This orientation was selected because it was 
initially thought that such an orientation would be convenient for the stri- 
ations model (Model No. 17). Later development indicated there was 
no need to have such an orientation. It can be easily changed; indeed, 
in the GRC version of the HAG module, the orientation has been made 


completely arbitrary. 


The azimuth (ALPHA) of the grid is set equal to the declination 
angle (DECANG) of the magnetic dipole field at the grid origin, obtained 
by calling Subroutine BFIELD. 


2.3 ESTABLISHMENT OF QUADRUPOLE COORDINATE SYSTEM 
FOR ARBITRARY QUADRUPOLE COLATITUDE AND 
QUADRUPOLE EAST LONGITUDE OF GRID ORIGIN 
Although we have described how the geographic azimuth of the 

quadrupole coordinate system has been established, it remains to de- 

scribe how the locations of the quadrupole-north and quadrupole-west 
poles are established. To do this we note that the geographic east 
longitude (GLO) and colatitude (GCO) of the grid origin are set by 

DRIVER reading Namelist SETUP1. The quadrupole east longitude 

(QLO) and colatitude (QCO) of the grid origin are also set, normally to 

90 deg but allowably in the range between 0 and 180 deg, by DRIVER 

reading Namelist SETUP1. These values of QLO and QCO, converted 

to radians and passed to Subroutine GRIDON through PARAMS Common, 
are used to establish the arc length (OQN) from the origin to the north 
pole in the quadrupole system, as illustrated in Fig. 2 and derived as 


follows. 
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Fig. 2. Geometry for Determining the Location of the Quadrupole 
North (QN) and Quadrupole West (QW) Poles. The quadrupole 
colatitude (QCO) and east longitude (QLO) of the quadrupole 
origin (O) are given. 
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For simplicity, assume 
0< QCO< 180 , 0 < QLO < 180 
By the law of cosines for angles, 
cos A - cos(QCO) cos(QLO) + sin(QCO) sin(QLO) cos 90 
= = cos(QCO) cos(QLO) 
By the law of sines, 


sin(OQN) = sin(QCO) sin 90 ’/sin A 


= sin(QCO)/(1 - cos“r)* 
or 
OQN = sin”! [sin(QCo)/(1 - cos” )] 
If 
QCO > 90° , 
set 


OQN = 7- OGN,, ; 


where ON, is the principal value of sin”! (sin OQN). Thus, the arc 
length OQN is determined and the quadrupole north pole (QN) is located 


with respect to the quadrupole origin O. 


We now must establish (a) the geographic east longitude (GLQN) 
of the quadrupole north pole and (b) the quadrupole east longitude (QLN) 
of the geographic north pole. These relations are needed in converting 


the coordinates of an arbitrary point from geographic to quadrupole 


coordinates and vice versa. 


SS ee 


By applying the law of cosines for sides to the triangle O-QN-N 


in Fig. 3a or Fig. 3b, we have 
cos(GCQN) = cos(GCO) cos(OQN) + sin(GCO) sin(OQN) cos|a@! , 
where 
0 < GCQN <= 180° , 
but we will need 
2 z 
sin(GCQN) = 1 - cos (GCQN)| 
By the law of sines, 
sin B = sin(OQN) sin! a@! /sin(GCQN) 
or 
é = sin®? |sin(OQN) sin| a| /sin(GCQN) 
From Fig. 3a (for a> 0) we see that 
GEQN = GLO+8B , a>0 


provided the angle § is not cut by the Greenwich (G) axis; if it is, then 
GLQN would exceed 27 and we would have to subtract 27 from GLOQN. 


From Fig. 3b (for a< 0) we see that 
GLQN = GLO-&8 , a<0 


provided the angle £ is not cut by the Greenwich (G) axis; if it is, then 
GLQN would be negative and we would have to add 27 to GLQN. Thus, 
GLQN is determined except for the effect of an uncertainty in the proper 


quadrant for p. 
To determine the quadrant for 8, we first consider the conse- 
quence of § equalling 90, i.e., by applying the law of cosines for sides 


to the triangle O-QN-N in Fig. 3a or Fig. 3b, we have 
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Fig. 3a. Geometry Used (for a> 0) in Deriving Geographic East Longitude 
(GLQN) of Quadrupole North Pole (QN) and Quadrupole East 
Longitude (QLN) of Geographic North Pole (N). Gis Greenwich. 
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Fig. 3b. Geometry Used (for a ~ 0) in Deriving Geographic East Longitude 
f (GLQN) of Quadrupole North Pole (QN) and Quadrupole East 
Longitude (QLN) of Geographic North Pole (N). G is Greenwich. 
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cos(OQN) = cos(GCO) cos(GCQN) + sin(GCO) sin(GCQN) cos 8 
cos(OQN), = cos(GCO) cos(GCQN) 


From Fig. 3a or 3b we see that if OQN > (OQN) , then B > 90°. How- 
aes the computer will return only the principal value 6, < 90°) of 
sin (sin 8). Thus, if B> 90°, we want B=7 - Ry Thus, the final 
formulas for GLQN are 


GLQN = GLO : e ‘ a20, Bs 90 


ceca alba ~, &z20, B> 90 


to which 27 is subtracted or added if GLQN is greater than 27 or less 


than zero, respectively. 


To obtain QLN we first apply the law of sines to triangle 
O-QN-N in Fig. 3a or 3b to get 


sin y = sin(GCO) sin! a! /sin(GCQN) 
y = sin’ [sin(GCo) sin! @| /sin(GCQN) 
From Fig. 39 (for a 0) we see that 
QLN = QLO = y + 27 for @>0Q 
and from Fig. 3b (for @ <0), 
QLN = QLO + y for a<0 


Thus, QLN is determined except for the effect of an uncertainty in the 


proper quadrant for y. 
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To determine the quadrant for y, we first consider the conse- 
quence of y equalling 90 , i.e., by applying the law of cosines for sides 


to the triangle O-QN-N in Fig. 3a or Fig. 3b, we have 


cos(GCO) cos(OQN) cos(GCQN) + sin(OQN) sin(GCQN) cos y 
cos(GCO),, = cos(OQN) cos(GCQN) 
From Fig. 3a or Fig. 3b, we see that if GCO > (GCO)p,, then y > 90. 
Thus, the final formulas for QLN are 
oa 
| | 
QLN = QLO - %o +27 for a>0O | 
if GCO< (GCO), 
= QLO + ee for «<0 | 


| QLN = QLO - (7 - Y») +27 for @>0O 
if GCO > (GCO), A 


= QLO + mots for a<o} 

Bae. ee en 
where %p is the principal value of sin’ tein y). 

2.4 CONVERSION OF GEOGRAPHIC COORDINATES OF A POINT 


TO QUADRUPOLE COORDINATES (SUBROUTINE GEOQUA) 


The conversion from geographic coordinates (GCX, GLX) of 
an arbitrary point to quadrupole coordinates (QCX, QLX), applied to an 
event-point X in Subroutine INDEX by a call to Subroutine GEOQUA, is 


derived as follows. 
From Fig. 4a (for € > 0) we see that 
QLXK = QLN - 27 +y , 


and from Fig. 4b (for ¢€ < 0) 
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Fig. 4a. Geometry Used (for € + 0) in Deriving Conversions Between 
Geographic (GCX, GLX) and Quadrupole (QCX, OLX) 
Coordinates of Point X. 


ee 


Fig. 4b. Geometry Used (for € <0) in Deriving Conversions Between 
Geographic (GCX, GLX) and Quadrupole (QCX, OLX) 
Coordinates of Point X. 
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QLX = QLN -y , e<0. 


From the law of cosines for sides applied to triangle X-QN-N in Fig. 4a 


or Fig. 4b, we have 
cos(XQN) = cos(GCX) cos(GCQN) + sin(GCX) sin(GCQN) cos 8 
but we need 
sin(XQN) = [1 - cos?(xQu)|* 


where 
B = GLQN - GLX , €>0 


= GLX - GLQN , oS | ae 


As can be seen from a polar-plot about the geographic north pole (N), 
8 cannot exceed 7. Thus, if from one of these expressions B is found 


to exceed 7, we must replace 8 by 27-8. From the law of sines, 
siny = sin(GCX) sin B/sin(XQN) 
or 


¥ = sin”? [sin(GCX) sin B/sin(XQN)] 


To determine the proper quadrant for y, we consider the con- 
sequence of y equalling 90; i.e., by applying the law of cosines for 
sides to the triangle X-QN-N in Fig. 4a or Fig. 4b, 


cos(GCX) = cos(XQN) cos(GCQN) + sin(XQN) sin(GCQN) cos y 


cos(GCX), = cos(XQN) cos(GCQN) . 


From Fig. 4a or Fig. 4b, we see that if GCX is greater than (GCX),,, 
then \|y| is greater than 90°. Thus, the formulas for QLX become 


ee RE ee 


QLX = QLN - 27 ee for «70 
if GCX < (GCX), 
= QLN - ¥ for «€- 0 | 
= QLN - 27 aula, for «> 0 


z GCX > (GCX), 


- QLN - (7 - yy) for € <0 | 


where Yy is the principal value of Sin tain y). 

It remains to establish formal criteria for determining the 
sign of angle «. From the polar plots in Fig. 5, one can establish the 
following criteria: 

For GLQN > GLX, 

€ 20 as GLQN - GLX ©. 17 
for GLQN < GLX, 

€=0 as GLX - GLQN — 7 


It is also useful to consider the diagram in Fig. 6 which shows values 
of € as the parameter 5 = GLQN - GLX varies. For 5 = 47, the 
angle y equals 0. For 6=0, y equals 7 if GCQN < GCX and 

y equals 0 if GCQN > GCX. 


Now we determine QCX. First, from the laws of cosines for 
sides applied to triangle X-QN-QW in Fig. 4a or 4b, we have 
cos(XQW) = cos 90 cos(XQN) + sin 90° sin(XQN) cos(QLX) 
= sin(XQN) cos(QLX) 


but we need 


2 


t 
sin(XQW) = [1 ~ cos xQw)| 
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a) e<0 
GLQN 4 GLX 


Ge— ~4-17-"f GLX 


QN 
«>0 e<Zy 
GLOQN < GLX 
Fig. 5. Polar Plots About Geographic North (N) Illustrating Dependence 


of the Sign of Angle « on the Relative Positions of Quadrupole 
North (QN) and the Point X. 
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Fig. 6. Illustration of the Dependence of the Sign of the 
Angle € onthe Parameter © ~ GLQN - GLX. 


From the law of sines, we have 
sin(QCX) = sin(KQN) sin(QLX)/sin(XQW) 


and 


QCX = in [sin(xQn) sin(QLX)/sin(KQW) 


To determine the proper quadrant for QCX, we replace Point O by 
Point X in Fig. 2 and note that QCX is less than 99° if XQN is less than 
90° and QCX is greater than 90’ if XQN is greater than 90°. 


2.5 CONVERSION OF QUADRUPOLE COORDINATES OF A POINT 
TO GEOGRAPHIC COORDINATES (SUBROUTINE QUAGEO) 
The conversion from quadrupole coordinates (QCX, QLX) of 
the arbitrary Point X to geographic coordinates (GCX, GLX), applied 
to the north and west sides of the columns in Subroutine GRIDON by a 
call to Subroutine QUAGEO although it is not currently used, is derived 


as follows. 


By applying the law of cosines for angles to triangle X-QN-QW 
in Fig. 4a or Fig. 4b, we have 
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cosA = - cos(QCX) cos(QLX) + sin(QCX) sin(QLX) cos 90” 


- cos(QCX) cos(QLX) 

but we will need 

sind = [1 - cos’) 
By the law of sines, we have 

sin(XQN) = sin(QCX) sin 90 /sin A 
= sin(QCX)/sin dA 
or 
XQN = sin”? ‘sin(QCX)/sin a] 


XQN is less than 90 if QCX is less than 90 and is more than 90° if 
QCX is more than 90 , as may be seen from Fig. 2. From Fig. 4a (for 


€>0), we see that 
Y QLX + 27-QLN , e>0 
and from Fig. 4b (for ¢€ <0), 
y = QLN - QLX , ex<U 4 


However, from the polar plots about quadrupole north in Fig. 7, we see 


that in general 
Y 'QLN - QLX! 


or 


y = 27 - |QLN - QLX| 


if |QLN - QLX| exceeds 7. By the law of cosines for sides applied to 
triangle X-QN-N in Fig. 4a and Fig. 4b, we have 


cos(GCX) = cos(XQN) cos(GCQN) + sin(XQN) sin(GCQN) cos y , 


or our desired expression for GCX, 


Fig. 7. 


en 


y = QLN - QLX| 


Y 


Qn - |QLN - QLX' 


Polar Plots About Quadrupole North (QN) Illustrating 
Dependence of Polar Angle y on Quadrupole East Lonmtudes 
of Geographic North (QLN) and Point X (QLX). 
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ca f 
GCX cos 1 |\cos(XQN) cos(GCQN) + sin(XQN) sin(GCQN) cos | : | 


We will also need 


i 2 


sin(GCX) i - cos“ (GCX) 


There is not a problem with respect to quadrant identification for GCX 
since the physical range of GCX, 0 to 7, is the same as the range of the 


principal values of cos” [cos(GCX) |. 


To determine GLX, we first apply the law of sines to triangle 
X-QN-N in Fig. 4a or Fig. 4b to get 


sin 8 = sin(XQN) sin y/sin(GCX) 
or 
Bp = ain’ |sin(XQN) sin y/sin(GCX) | 


To determine the proper quadrant for 8, we consider the consequence 
of 6 equalling 90 , i.e., by applying the law of cosines for sides to 


triangle X-QN-N, we have 
cos(XQN) = cos(GCX) cos(GCQN) + sin(GCX) sin(GCQN) cos 8 
cos(XQN), = cos(GCX) cos(GCQN) 
Thus, § is less than 90° if XQN is less than (XQN), and 6 is greater 
than 90° if XQN is greater than (XQN),. 


Finally, from Fig. 4a (for € > 0) and Fig. 4b (for € < 0) we 


see that 


GLX = GLQN *B, for € 20, XQN < (XQN), 


GLX = GLQN ¢ (7 - B,) for € <0, XQN 2 (XQN) 2 , 
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where A, is the principal value of Sin etn B). 


It remains to establish formal criteria for determining the sign 
of the angle €. From the polar plots in Fig. 8, one can establish the 
following criteria: 

For QLN > QLX, 
€=0 as QLN - QLX 7 ; 
for QLN < QLX, 
€=0@ as QLX - QLN -.7 
It is also useful to consider the diagram in Fig. 9 which shows values 
of € as the parameter  ~ QLN - QLX varies. For = +7, the angle 


B equals 0. For yp =0, 8 equals 7 if QCN < QCX and £8 equals 0 if 
QCN > QCX. 


2.6 SPECIAL TRANSFORMATION PROBLEMS 


If the quadrupole north is in close proximity to geographic 
north, the transformation of the colatitude of Point X from one coordi- 
nate system to the other is calculated in the usual manner. However, 
there may be an angular displacement in the longitude of Point X between 
the two coordinate systems. This angular displacement is the difference 
between GLO (geographic longitude of the origin) and QLO (quadrupole 
longitude of the origin). Thus, in transforming from the geographic to 


the quadrupole system, the quadrupole longitude of Point X is given by 
QLX = GLX - (GLO - QLO) 


in transforming from the quadrupole to the geographic system, the geo- 


graphic longitude of Point X is given by 


GLX = QLX - (QLO - GLO) 
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Fig. 8. Polar Plots About Quadrupole North (QN) Illustrating 
Dependence of the Sign of the Angle € on the Relative 
Positions of Geographie North (N) and the Point X. 
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Fig. 9. Illustration of the Dependence of the Sign of the 
Angle € onthe Parameter pp = QLN - QLX. 


In the general situation, where quadrupole north and geographic 
north are separated, there is one caution to be noted. If the Point X 
lies on the great circle passing through both quadrupole north and geo- 
graphic north, the point may not be transformed correctly from one 


system to the other without a modification of the code. 


2.1 QUADRUPOLE COORDINATES OF TOP-LEFT CORNER 
OF GRID 


The quadrupole east longitude (QLMIN) and colatitude (QCMIN) 


of the top-left corner of the grid (see Fig. 10) are computed from 
QLMIN = QLEBOT + 0.5 x DELY 
QCMIN = QCSBOT + 0.5 * DELX 


where 
QLEBOT = QLO - NYQW * DELY 


QSCBOT = QCO - NXQN * DELX 


and NYQW and NXQN are the number of columns to quadrupole west and 
north, respectively, of the grid origin and DELX and DELY are the 
angular widths of the columns in the quadrupole X-direction and 


Y-direction, respectively. The quantities NXQN, NYQW, DELX, and 
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Quadrupole Grid. Schematic plan view of the quadrupole grid 
shows location of the top-left (north-west) corner with 
respect to grid origin O. Cell boundaries are not squares 
(or even rectangles) as shown for simplicity but are great 
circles. 


DELY are set by DRIVER reading Namelist SETUP1 and passed to 
GRIDON through PARAMS Common. 


2.8 COORDINATES OF THE NORTH AND WEST SIDES OF THE 
COLUMNS 


2.8.1 Quadrupole Coordinates 


The quadrupole coordinates, east longitude (QL(J)) and colati-~ 
tude (QC(I)), of the west and north sides of Column (I, J) are computed 


from 


1, NCELX 


QC(I) = QCMIN + DELX x (I- 1) , I 


QL(J) = QLMIN + DELY x (J - 1) , J =1, NCELY 


where I and J are the column indices in the quadrupole X-direction and 
Y-direction, respectively, and NCELX and NCELY are the total number 
of columns in the quadrupole X- and Y-directions, set in DRIVER and 
passed to GRIDON through PRPREG Common. The Sines and cosines 
of both QC(I) and QL(J) (SNQC(I), CSQC(I), SNQL(J), and CSQL(J)) are 


also computed. 


2.8.2 Geographic Coordinates 
DO-loop 32 in GRIDON, developed as a possibly-needed capa- 


bility to convert column coordinates from the quadrupole coordinate 
system (QC(I), QL(J)) to the geographic coordinate system (GC(I, J), 
GL(I, J}} is currently bypassed. The transformation is derived in Sec- 
tion 2. 5. 


40 


eee 


2.9 QUADRUPOLE COORDINATES OF SOUTH AND EAST SIDES 
OF THE COLUMNS 


The quadrupole coordinates, east longitude (QLE(J)) and colati- 
tude (QCS(I)), of the east and south sides of the Column (I, J) are com- 


puted from 


QCS(I) = QCS(I - 1) + 2 x [QC(I) - QCS(I - 1)] I=1, NCELX 


QLE(J) = QLE(J - 1) + 2 [QL(J) - QLE(J - 1)] 


oy 
| 


= 1, NCELY 
where 

QCS(0) = QCSBOT 

QLE(0) = QLEBOT 


The sines and cosines of both QCS(I) and QLE(J) (SNQCS(I), CSQCS(I), 
SNQLE(J), and CSQLE(J)) are also computed. 


2.10 SOLID ANGLE OF COLUMN IN QUADRUPOLE 
COORDINATE SYSTEM 


2.10.1 Background 


Whereas the actual solid angle subtended by a column of cells 
is not needed for computing the hydrodynamics per se, there are at 


least two requirements for such information: 


a. The actual mass ina cell is needed tc compute the specific 
energy (erg/g) deposited by heavy particles from the loss- 
cone and ion-leak sources [Vol. 5, HS-75b|, and 


b. The mass of the event cell is required to comp*re the total 
energy put into the event cell with the total erergy of the 
event. 


ee oe 


2.10.2 Derivation 
The element of area on a sphere formed by the quadrupole 
planes (@, @ + da, B, B + dB), as shown in Fig. 11, is a parallelogram 
with sides rd@ and rdg and with acute angle ©, given by 
Os Mik 
where A is given by the law of cosines as 
cos A - cos &@ cos Bf + sin @ sin B cos 90 


= =- cos @ cos Bp 


Hence, 


cos 59 = cos(m - A) = - cosA = cos a@ cos B 


4 4 
é é 


sin 0 = (fh - cos“6) =(1- cosa cos?) 


Thus, the solid angle for quadrupole angles OD Os By Bo is 


a B = 
Q(a,, Oo) By Bo) = f aff (1 - eos “cr cos"p)” da dp 
sc | 
Since the first integration leads to elliptic integrals (of the second kind), 
we have resorted to the following approximation, the validity of which 
can be checked by numerical integrations, if necessary: Since 


a~B~1/2, we know cos“a cos“B < 1; hence 
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Geometry Used in Computing Solid Angle of Column in 
Quadrupole Coordinates. 


é. 8 $ 
81(4) 5455 By» By) i 2 , 2 ( - 5 cosa cos“8 da dg 
1 1 


= (a, - a4) (By - B,) - E (a, - a, + sin A cos om) 


- sin a, cos a,) (By * By + sin By cos By - sin B, cos B,) 


i 


1 1k ae : 
(a, - a4) (By - B,) = 5 la S a, + 3 (sin 20, ~ sin 204) 


x lbp = By ¥ 5 (sin 285 - sin 28.) 


The solid angle of each column is computed and stored as a 
one-dimensional array, STR(NN), in PARAMS Common. 


2.11 CELL BOUNDARIES (FOR SAIHYD) 


The bottom boundary of the grid (HABOT), set by DRIVER 
reading Namelist SETUP1, is passed to GRIDON through PARAMS Com- 


mon. The (initial) bottom boundary of Cell K is computed from 


HB(K) = HBO(K) = HBO(K - 1) + DELHA(K - 1) , K = 2, KMAX 
where 

HBO(1) = HABOT 

KMAX = NCELZ +1 


DELHA(K), the initial height of Cell K, is set by DRIVER reading Name- 
list SETUP1 and passed to GRIDON through PARAMS Common. NCELZ, 
the total number of cells ina column, is set by DRIVER reading Name- 
list SETUP1 and passed to GRIDON through PRPREG Common. 


The geocentric radius of the bottom boundary of Cell K is com- 
puted from 
R(K) = RE + HB(K) 


2.12 STERADIANAL VOLUME OF CELL (FOR SAIHYD) 


The steradianal volume of Cell (K - 1) is computed from 


VOL(K - 1) = {[R(K) P - ([R(K - 1)P}/3 


2.13 CELL-CENTER ALTITUDES (FOR SAIHYD) 


The altitude of each cell center, including a fictitious cell- 


center above the top boundary, is computed from 
HC(K - 1) = 0.5 x [HB(K - 1) + HB(K)] , K = 2, KMAX 


HC(KMAX) = 0.5 x [3 x HB(KMAX) - HB(NCELZ)| . 


3. SAITHYD (DIFFERENCE-EQUATION FORM OF 
LAGRANGIAN HYDRODYNAMICS) 


3.1 INITIALIZATION OF LAGRANGIAN CELLS IN A COLUMN 
UNDER AMBIENT CONDITIONS 


The initialization of the Lagrangian cells in a column is started 
by a call from Subroutine GRIDON to Subroutine ATMOSG. 
ook Velocity 


The velocity of the bottom boundary of each Cell K, VEL(K), 


is set to zero: 


VEL(K) = 0 , K=1, KMAX 


3.1.2 Viscous Pressure 
The viscous pressure of each Cell K, Q(K), is set to zero: 


Q(K) = 0 , K=1, NCELZ 


Orbea Mass 


The mass per steradian of Cell K is found by a call to Subrou- 
tine SRMASS which uses the Weddle Rule to integrate the function oR” 
between the lower and upper geocentric radii R(K) and R(K+1) for 
Cell K, i.e., 


7 
CELM(K) = (0.3/6)[R(K + 1) - R(K)] > p, R? w, 
i=l 


where W, = 1, 5, 1, 6, 1, 5, 1. The mass density p; is found at each 


of the seven evenly-spaced geocentric radii R. by calls from Subroutine 
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SRMASS to Subroutine ATMOSU. The mass of fictitious Cell KMAX, 
CELM(KMAX), is arbitrarily set equal to CELM(NCELZ). 


The so-called boundary mass is computed as the mean of the 


two adjacent cell masses, i.e., 


BDYM(K + 1) = 0.5 x [CELM(K) + CELM(K + 1)] , K = 1, NCELZ 


3.1.4 Pressure , 


The pressure in Cell K is computed from the hydrostatic form 
of the momentum equation so as to stabilize the column of cells under 
ambient conditions. (This procedure is patterned after that in the MRC 
HARC code.) We start with the ambient pressure (PAMB) at the 
fictitious-cell center above the top boundary, obtained by a call to 


ATMOSU at altitude HC(KMAX), and work downward, using the expression 


P(K - 1) = P(K) + G[R(K)]* BDYM(K)/[R(K) |? , 


where G|R(K)]| is the acceleration due to gravity at geocentric radius 
R(K), 


(2 


G|R(K)| = GRAVZ ™ }RE/[RC(K)|\", 


and GRAVZ, the value of G at the Earth's surface, is set in Block Data 
BLOCKH and passed to ATMOSG through CNSTNT Common. 


3.1.5 Density 


The mass density of Cell K, RHO(K), is computed from 


RHO(K) = CELM(K)/VOL(K) 


3.1.6 Internal Energy 


The specific internal energy of Cell K, EPS(K), is computed 


from the equation of state, 
EPS(K) = P(K)/|RHO(K) = GAM1| , 


where GAM1 = y - 1 and is set to be 0.5 in Block Data BLOCKH and 
passed to ATMOSG through CNSTNT Common. 


nl a Species Densities and Other Chemistry Quantities 


Species densities and other chemistry quantities are not re- 
quired for the hydrodynamics but they are required for the chemistry. 
In general, these quantities are determined in Subroutine ATMOSG at 
the beginning of a loop over the cells in a column by three calls to Sub- 
routines ATMOSU(2, HCKM), IONOSU(2, HCKM), and SPCMIN(2, HCKM), 
where HCKM is the altitude (in kilometers) of the cell center. We list 
in Table 1 the chemistry quantities and how they are initialized; this 


table applies to NRLHYD as well as SAIHYD. 


3.2 CONTROL ROUTINE FOR ADVANCING THE HYDRODYNAMICS 
(SUBROUTINE HYDROG) 


Subroutine HYDROG, the routine controlling the subroutines 
doing the actual advancement of the hydrodynamics, performs the fol- 
lowing functions. 

a. Establishes the proper headings for the two sets of time- 

slot data by testing on the value of the parameter (IDREG) 
for the number of hydro and chemistry updates. IDREG, 


initially set and updated in DRIVER, is passed through 
PRPREG Common. 


b. Controls the looping over all the columns in the grid 
while the hydrodynamics is being done. 


Table 1. Initialization of Chemistry Quantities for 
SAIHYD and NRLHYD. 


Index M 


in Fortran Initial 
BUF(L+M) Quantity Variable Value Basis 

TX Formula 
ENE EFE x FP 
EN4S SNI(7) x F 
EN2D 1.0 
SNI(3) SNI(3) * F 
SNI(1) SNI(1) x F 
SNI(2) SNI(2) x F 

{NO | ENO SNI(8) * F 

[N° | ENP 1.0 

[o"| OP EFOP x F 


TNovib TN2VIB TX 


[CO, | SNI(6) SNI(6) * F 
Q, mb EFQ EFQ 
|He | SNI(5) SNI(5) x F 


[M "| | AMOLP | Max(1, e-N’-O°) 


4similar to formula in Subroutine IONOSU ° Vol. 14a) but with account of 
normalized species. 

by, known as F ACT in code, is the factor required to normalize the species densi- 
ties so that the sum of their mass densities will equal the mass density computed 
for the Lagrangian cell (see Section 3. 1. 5) for SAIHYD or for the Lagrangian point 
(see Section 4.1.4) for NRLHYD. 

Subroutine IONOSU through IONOUP Common. 

Subroutine SPCMIN through ATMOUP Common. 

€Arbitrary. 

fSubroutine ATMOSU through ATMOUP Common. 

&STN2VIB = TX in ambient ionosphere model. 
Not a stored quantity for either SAIHYD or NRLHYD but is computed for 
information. 

‘Charge conservation. 
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c. Effects the transfer of the NCOL cell quantities for the 
column of interest from large core memory (LCM) to small 
core memory (SCM). This transfer is done by a call to 
Subroutine ECRD(BUF, IWI, NCOL), where IWI = (NN-1) 

« NCOL + NTIME1, so that NCOL quantities starting at 
index TWI in LCM are transferred to the BUF array in 
SCM. 


d. Initializes the variables in CELLLS Common and HEICiIT 
Common before calling Subroutine HYDRO1. These vari- 
ables for Cell K (K = 1, NCELZ) are: 


HC(K) = BUF(LOC + 1) Cell-center altitude 
RHO(K) = BUF(LOC + 2) Mass density 
VEL(K + 1) = BUF(LOC + 6) Upper-boundary velocity 
P(K) = BUF(LOC + 15) Pressure 
Q(K) = BUF(LOC + 16) Artificial viscous 
pressure 
HBO(K + 1) = BUF(LOC + 21) Top-boundary altitude in 


ambient atmosphere 
R(K + 1) = R(K) + 2 * |RE + HC(K)- R(K)| Geocentric radius of top 
boundary 
BUF(LOC + 15) 


EPS(K) = BUFILOC . 2) = GAMI Specific internal energy 
VOL(K) = || R(K + vp ~ |R(K) P1/3 Steradianal volume 
CELM(K) = RHO(K) * VOL(K) Cell mass 
CELM(KMAX) = CELM(NCELZ) Mass of fictitious top cell 
PKMAXZ = PKMAX(NN) Pressure in fictitious 
top cell 


BDYM(K + 1) = 0.5*|CELM(K)+CELM(K+1)| Upper boundary mass 


e. Calls Subroutine HYDRO1(TIME, TSTOP, IEXPLD) to ad- 
vance the hydrodynamics in a column from the start time 
(TIME) to the stop time (TSTOP). A value of one for 
IEXPLD, set in DRIVER and passed to HYDROG as a calling 
argument, denotes that an event has just occurred and a 
value of zero denotes that an event has not just occurred. 

A value of one for this flag is used in Subroutine DELTIM 
to set a small value for the first time step. This flag was 
used in development and may no longer be necessary. 


eA 


—~ 


‘1. Checks whether or not Column NN needs to be rezoned, 


based on the value of the geocentric radius of the topmost 
boundary relative to the criterion radius RMAX, and sets 
the rezone flag IREZ(NN) if required. IREZ(NN) is incre- 
mented by one if rezoning is required for Column NN, thus 
making its value odd. After rezone is performed in Sub- 
routine REZONE, based on IRE Z(NN) being odd, IREZ(NN) 
is then incremented by one in Subroutine REZONE, thus 
making its value even. The number of rezones already 
performed on Column NN can be found by dividing IREZ(NN) 
by two. 


.2. Prints a message, if a rezone is required, stating the re- 


zone number and the associated time. 


g. Outputs the following hydrodynamic quantities for Column 
NN at time TSTOP: 


Cell index 

Cell-center altitude 

Cell mass density 

Cell-center velocity 

Cell specifie internal energy 

6. Cell pressure 

7. Cell artificial viscous pressure 


Ton Fe wh 


h. Prints the quantities ETZ, ET, and FDELE from the energy 
check, computed in Subroutine ENECHK after being called 
by HYDRO1 and passed through ENERGY Common. 


i. Updates the BUF array. (Compare Step d above.) 


j. Effects the transfer of the NCOL cell quantities for the 
column of interest from SCM to LCM by a call to Entrv 
Poiat ECWR(BUF, IWI, NCOL) in Subroutine ECRD. 


k. Sets and prints the values of the altitudes and velocities of 
the time-dependent event-points by calling Subroutine 
TIMVAR. 


mu 


3.9 HYDRODYNAMICS EQUATIONS AND TIME STEP 

(SUBROUTINES HYDRO1 AND DELTIM) 

In SAIHYD, the basic hydrodynamics is computed in Subroutine 
HYDRO1, which advances the hydrodynamics in a geocentric steradianal 
column of cells by solving the one-dimensional Lagrangian hydrodynam- 
ics equations in a geocentric spherical coordinate system. The differ- 
ence equations and their sequence of solution are patterned generally 
after those used by D. Sappenfield in his version of RADFLO |MRC-R- 
67 (generally unavailable) and computer listing kindly provided on 
18 April 1974 to D. A. Hamlin], originally developed by J. Zinn |Zi-73a| 
of Los Alamos Scientific Laboratory. We elected this procedure in 
preference to using the MRC HARC-Code procedure which had undergone 


an obscure evolution from Sappenfield's. 


The changes we have made in Sappenfield's procedures include 
the following: 


a. The equations are solved for flow in one steradian instead 
of 47 steradians. 


b. The acceleration due to gravity is added. 


c. The time step (DT) is determined, as in the MRC HARC- 
code, by solving the White stability condition for DT, 
instead of computing the White number and changing DT 
according to some prescription. 


d. The total artificial viscous pressure is currently taken to 
be a sum of the linear and quadratic viscous pressures, 
each weighted by multipliers EPS1 and EPS2 (read by 
DRIVER in Namelist SETUP1 and passed to HYDRO1 
through HEIGHT Common), whereas Sappenfield uses only 
the linear form. We are currently using a value of 0.5 for 
both EPS1 and EPS2, but this topic needs more study. 
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In Sappenfield's equations, the velocity is not at time ae as 


« 


in many or most explicit (as opposed to implicit) Lagrangian hydro 
codes, but is at time ts just aS are the radius, energy, etc. In com- 
puting the new velocity, Sappenfield uses an acceleration which is so 
chosen that the sum of the kinetic and internal energies is conserved to 


second order in the time step, aS we now demonstrate. 


For simplicity, we consider two cells in a one-dimensional 
cartesian system, illustrated in Fig. 12. Assume the boundaries are at 
rest, i.e., bs ae 0. At the beginning and end of a time step 


(At =t eit 


) the kinetic energies, KE(t_) and KK(t ), are 
n+1 n n n+1 


m + Mm 
; a ae BV ie 
KE(t,) = 5 (+) [Vv 


1 at + te ; D) 
y= Ee trot) + ae 


KE(t oy 5 


where (my : M9) ‘2 is the boundary mass and @ is the acceleration. 


The increase in kinetic energy, AKE, is 


AKE = KE(t. .)}- KE(t,) 


n+l 


LS) el 


m, +m 
1 2 2 
(1? [2 Volt.) a At + (@ At) 


1 sa + Ths 1 
LY z— [vott,) 50 at] Qa At 


The change in internal energy, OIE, is 


AIE - (Py AV, +p, AV 


ee Reg, V. 


ty Vo 
ear, oN 
| 1 


Fig. 12. Two Cells ina One-Dimensional Plane-Geometry 
Lagrangian Hydro Code. 


where the change in volume is 


Thus the change in internal energy 


iS 


AIE AV , (Po - Py) 


and 


oe 
AKE 


If we want AIE to equal AKE, 


ye 


1) 


a At| At(p, - p 


must use 


{ ) 


an acceleration @ given by 


ewe 


me. 


which is the form of the acceleration used by Sappenfield. 


The time step is normally determined by calling Subroutine 
DELTIM. Before ‘deriving’ the expression given in DELTIM for DT, 
we first recall the Courant condition for determining the stability of the 
difference equations of hydrodynamics. The Courant condition |RM-67, 
p. 298|, 

c 
DELR/DT 
where c is the sound speed and DELR is the cell width, states |Hi-65a | 
that "the time step must be small enough so that a disturbance traveling 
at local sound speed cannot travel a full zone width, otherwise hydro- 


dynamic disturbances can escape consideration." 


George N. White, Jr., of LASL has long used a generalization 
of the Courant condition for determining the stability of the difference 
equations of hydrodynamics. His method was initially derived empiri- 
cally for use in early Los Alamos codes and later derived in a memoran- 
dum, entitled "An Analysis of Lagrangian One-Dimensional Shock 
Hydrodynamic Difference Equation Stability by Application of Routh's 


Criterion," kindly provided by C. G. Davis of LASL. 


This memorandum had been provided to R. D. Richtmyer, for 
Richtmyer and Morton |RM-67, addendum on p. 350] state that G. N. 
White communicated an extensively-used stability condition for the vis- 


cosity method: 


- ¥ o\¥ry - Veal 
AX/ At n 
ae 
fu: 
a 


where we have corrected a misprint in RM-67 by writing the first term 
oF ‘ : 2 
as a squared quantity instead of as a first-power quantity. Here, a is 
the dimensionless coefficient in the quadratic form of the artificial 
viscosity. 
To show the expression we use for the time step, we start 


with the expression 


2 \ 
DT DDTSRV * DT 
WN (oral cs2.4 Dees (1) 


DELR * (WS) 
where 
WN = the White number, set to be 0.6 in Subroutine DELTIM. 
Di LR = RTOP - RBOT 
WS = (RTOP + RBOT)/2 


CS2 = alleged ''sound-speed squared" 
= 4.166 P/p 


DDTSRV = rate of change in steradianal volume. 


DDTSRV * DT = change in steradianal volume. 

DELR » (ws)? = approximate steradianal volume. 

Equation (1) is similar to Eq. (9) in G. N. White's memorandum; how- 
ever, it differs in that (a) the details of computing the volume are 
changed, just as RM-67 altered them, (b) the sound-speed squared has 


the mysterious factor 4.166, unrecognized by Sappenfield, Zinn, or 


MRC personnel, and (c) White's a? has been set equal to 1. Equation (1) 
is very close to Sappenfield's expression for the White condition; it dif- 
fers in that (a) Sappenfield's first term is only one-half as large and 


(b) the details of the volume computation differ. 


56 


ete 


oe 


Now, rewrite Eq. (1) as 


1 = (2 y Cs2__4 DDTSRV x DT 
DELR) WN * WN pore x ws)? 


and define 


2 
| 


B = 2 DELR X DDTSRV/|CS2 » (WS) 


C = (DELR)2 WN/CS2 


so that 
Gee (DT)” , DT x 2B 
ais 5 C 
or 
(DT)? + 2Bx DT-C =0 
or 


which is the expression in Subroutine DELTIM for determining the time 


step. 


Having the time step, we are now ready to advance the hydro 
equations in time, done by looping over the cells in a column in two 
phases. In Phase 1, we compute the total pressure in each of the cells 
at time ‘. by adding the real pressure and artificial viscous pressure 
(from the previous Phase 2). The real pressure is computed by using 
a constant-gamma equation of state, with the mass density and specific 
internal energy from the previous Phase 2. Thus, the total pressure 


in Cell K is 


PTOT(K) = P(K) + Q(K) K= 1, NCELZ 
57 
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ee ee 


Our assignment of the boundary-condition pressure at the fictitious 

cell center above the top boundary is an ad hoc procedure that may need 
improvement. Currently, this constant pressure (prior to a rezone) is 
taken to be the ambient pressure at the initial altitude of the fictitious 


cell-center above the top boundary: 
PTOT(KMAX) = PKMAXZ , 


where PKMAXZ is determined in Subroutine ATMOSG and passed to 
HYDROt1 through HEIGHT Common. 


In Phase 2, the difference equations, applied in sequential order 
to each cell before moving on to the next cell, yield values for the fol- 


lowing quantities: 


a. Mean acceleration of boundary during DT. 


b. Radius of boundary at time tet’ 


c. Geocentric radius (and altitude) of cell boundary at time 


Le 


d. Thickness of cell at time tet’ 


e. Velocity of cell boundary at time tad’ 


~ 
. 


Steradianal volume of cell at time ty 1 


Change in steradianal volume of cell during DT. 


a 


h. Specific internal energy of cell at time tel 


_ 
. 


Mass density of cell at time t : 
n+1 


j. Artificial viscous pressure of cell at time to. 1 


We let K = 1 denote the lowest boundary and lowest cell, and 


K1 = K- 1. Then the equations are: 


t 
| Acceleration: 
“cyP — _§[R(K) 2 . DT x R(K) VEL(K)| 
ACC = |PTOT(K1) - PTOT(K)|> BDYM(K ‘7 G| R(K) | 
58 
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Re eR 


where 


G(A) = GRAVZ ™ (RE/A)? 


R(K) = R(K) + [VEL(K) + 0.5 x ACC x DT| x DT 
Altitude: 
Boundary: 
HB(K) = R(K) - RE 
Cell Center: 
HC(K1) = 0.5 * |HB(K) + HB(K1)| 
Cell Thickness: 
DR = R(K) - R(K1) 
Velocity: 
VEL(K) = VEL(K) + ACC x DT 
Volume of Cell Kt: 
VOLUME = /[R(K)P - [R(K1) [t/3 
| Volume Change of Cell Ki in DT: 
| DVOL = VOLUME - VOL(K1) 
Store Volume of Cell K1: 
VOL(K1) = VOLUME 
Internal Energy of Cell K1: 


ITOT(K1) * 
EPS(K1) = EPS(K1) - PPORCSE VOL 


Density: 
CELM(K1 
REOUED) TOLKTY 


Artificial Viscous Pressure: 


| EPS1* QL+ EPS2xQ2, DVOL<0 
Q(K1) 
lo .  DVOLz0 


where the linear and quadratic viscous pressures are 


Qi = - [1.67 * PTOT(K1) x RHO(K1)* x 0.4 * DDDV 


Q2 = 0.4 * RHO(K1) x (DDDV)” 


DDDV = DVOL x DR/|DT X VOL(K1)| 


Real Pressure: 


P(K1) = GAM1 * RHO(K1) * EPS(K1) 


3.4 ENERGY CHECK (SUBROUTINE ENECHK) 


We have provided an energy check for the SAIHYD, to be used 
at least during the development stages of the code. Such a check is 
made in Subroutine ENECHK. This routine is called to perform one of 
two types of calculations: (1) To compute the total energy (ETZ) ina 
column, consisting of internal and kinetic energy, regarded as the refer- 
ence energy for comparison purposes and (2) To compute the current 
total energy (ET) ina column, consisting of not only the current internal 
and kinetic energies but also the total work done by the column since the 
reference energy (ETZ) was last computed and updated (in Subroutine 


CHEMG) by the return of energy from the chemistry routines. The total 
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work done consists of (a) work against gravity and (b) work by the top 
cell against the downward-acting backpressure. After ET is computed, 
the fractional change in the total energy (FDELE) is computed, 
FDELE = (ET - ETZ)/ETZ 
The formulas used are, for K= 1, NCELZ: 
Internal Energy: 
EI = >, [EPS(K) x CELM(K) | 
K 
Kinetic Energy: 
EK = 0.5 > =. BDYM(K+1) » [VEL(K+1)|? 
K 


Work Against Gravity: 


aE ee , A 1 i 
Bis eee Bas STE ” RE + 0.5[(HBO(K) + HBO(K+1) RE + HC(R)| 


Work Against Backpressure: 


EW = PKMAXZ * |(RCKMX)® - (RCKMXZ)°|/3 


where 
EPS(K) = specifie internal energy (erg/g) 
CELM(K) = cell steradianal mass (erg/sr) 
BDYM(K) = boundary mass (g/sr) 
VEL(K) = boundary velocity (cm/sec) 


2 


GZREZ = ¢g_ (RE) (em*/see*) 


19] 


PKMAXZ = constant backpressure (dyne/em?) 
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RCKMX = RE + HC(KMAX) 
RCKMXZ = RE + HCKMXZ 


HCKMXZ = 0.5 [3 HBO(KMAX) - HBO(NCELZ)| 


A zero (or unit) value for an index flag (IENCHK), read by DRIVER 
with Namelist SETUP1 and passed through LINK Common, is used to by- 
pass (or permit) the use of Subroutine ENECHK. Within Subroutine 
ENECHK, a zero value for the hydro-cycle counter in Column 1 (KCYC(1)) 


allows ETZ to be computed and a nonzero value allows ET to be computed. 


Subroutine ENECHK computes ETZ upon calls from (a) DRIVER 
if ambient conditions are being evaluated or if energ, las been deposited 
for an event and (b) Subroutine REZONE after the column has been re- 
zoned. Subroutine ENECHK computes ET upon a call from Subroutine 


HYDRO1 at the end of a hydro update interval. 
3.5 REZONE (SUBROUTINE REZONE) 
de Dek Purpose 


To provide for maintaining a reasonable cell resolution in the 
(more important) lower portion of the HAG as a column expands to very 


high altitude, we have developed Subroutine RE ZONE. 


KG a General Features of Subroutine REZONE 


Subroutine REZONE resets the boundaries and properties of the 
cells for those columns whose geocentric radius of the uppermost bound- 


ary exceeds a given geocentric radius RMAX at the hydro-advance time 


TSTOP. The number of cells in a rezoned column is kept the same 
(NCELZ). The new cell thicknesses (DELHA2(K)) are read by DRIVER 
in Namelist SETUP1 and passed through REZONN Common. If 


DELHA2(1) = 0., the rezoned thicknesses are set equal to the original 
cell thicknesses (DELHA(K)) set in DRIVER by a Data Statement. Ta- 


ble 2 gives the currently-used values for the initial and rezoned cell 


thicknesses and altitudes. 


In our original version of Subroutine REZONE, we used the 
principles of conservation of mass, momentum, and energy to derive 
the rezoned quantities. This version usually worked well, but we some- 
times had (1) a negative internal energy, presumably due to the coarse- 
ness of the cell structure, and (2) one or more new boundaries with zero 
velocity, if the lowest old cell expanded so much that its center became 
higher than the upper boundary of one or more new boundary-cells. 
Failing to find any coding errors, we reluctantly gave up the principles 
of momentum and energy conservation and resorted to linear interpola- 
tion to find the new velocities and logarithmic interpolation to find the 
new pressures. However, we have retained conservation of mass within 
the boundaries of each new cell and a similar conservation of the total 


electron thermal and excitation energy. 


Subroutine REZONE is called from DRIVER after all columns 
have been advanced first for hydrodynamics and then for chemistry. In 
looping over the columns, Subroutine REZONE rezones only those col- 
umns for which a rezone flag (IREZ(NN)) has been set (based on the re- 
zone criterion) in Subroutine HYDROG immediately after Subroutine 


HYDRO1 completes the hydrodynamics for Column NN. 


If the rezone flag has been set for Column NN, then Subroutine 
ECRD is called to obtain the column quantities from Time-Slot 2. For 
development purposes we print the status of the Lagrangian-cell quanti- 


ties, and also the molecular-ion density and heavy-particle temperature, 


before rezone. After the new cell quantities are determined, Entry Point 
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ECWR in Subroutine ECRD is called to store the cell quantities in Time- 
Slot 2. Again, for development purposes, we print the rezoned cell 
properties, including the molecular-ion density and heavy-particle 


temperature. 


5 Fats A: Altitude in Ambient Atmosphere Corresponding to New 

Cell Locations 

One of the needed properties of a cell is the altitude (HBO(K)) 
in the ambient atmosphere corresponding to the lower boundary of the 
cell. Such information is used by Subroutine ENECHK in computing the 
work done against gravity since the start of the problem or since the 
last rezone. Subroutine REZONE uses linear interpolation between the 
appropriate pairs of the old values of these ambient-atmosphere altitudes 


to obtain the new values. 


3.5.4 Conservation of Mass, Species, and lon-Production Rate 


To facilitate implementation of the conservation laws, we de- 
veloped Subroutine FRCVOL, called from Subroutine REZONE, to com- 
pute the (absolute) fractional volumes (DVOL(NI)) of a new cell sliced 
by old-cell boundaries. In this way, e.g., the total mass of a new cell 
is equal to the sum of the products of the mass density and fractional 
volume for the old cells overlapping or contained within a new cell. The 
density in the new cell is obtained by dividing its total mass by its total 


volume. 


A similar procedure is followed for each density-dependent 
quantity, e.g., the species number densities and the effective ambient 


ion production rate. 


NS Pare abo. Conservation of Sum of Eleetron Thermal and Excitation 


Energies 


To compute the new value for the temperature taken to be the 
common value for the electrons, N, vibrational level ind the o('p)- 
to-O(°P) population ratio, we conserve the sum of the electron thermal 
energy and the excitation energy for op) and Ny vibration In com- 
puting this excitation energy we use the same parameters as used in th 
GET-~-developed Subroutine TEXK |KJ-74, KJ-74b|. After the amount 
of energy per unit volume available for partitioning among the three 
modes is known, as well as the densities of N, molecules, O atoms, ay d 


electrons, Subroutine TEXK is called from REZONE to return the te) 


perature common to the three modes. 


3.5.6 Velocity Interpolation 


The velocity of a new cell boundary is found by linearly inter- 


polating between the velocities of the two bracketing old boundaries. 


is Dai k Pressure Interpolation 
The pressure of a new cell, assumed to be specified at the 

center of the cell, is found by linearly interpolating between the loga 

rithms of the pressures at the two bracketing old-cell centers. In this 


procedure, the pressure (PKMAX(NN)) in the fictitious « ell at the top of 


the column is determined in the same way. 


3.5.8 Artificial Viscous Pressure 
In HYDRO1 we compute the artificial viscosity from 
Q(K1) = EPS1 * QI + EPS2 x Q2 (2) 


where 


Q1 - SQRT(1. 67 * PTOT(K1) * RHO(K1)) * WB © DDDV (3) 
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Q2 WB * RHO(K1) * DDDV * DDDV 
PTOT(K) P(K) + Q(K) 
WB 0.4 
DDDV DVOL * DR/(DT * VOL(K1)) 
DVOL 


(Volume) - (Volume) 
new 


old 


There are two difficulties in computing Q upon a rezone: 


1. How to compute DDDV ? 


2. How to handle PTOT (since it depends on Q)? 
These difficulties have been overcome: 


1. To handle DDDV, we replace DVOL/DT in (7) by 


DVD: R2 XU. - R2 U ; 
i i i i-1l i-1l 
so that 
DDDV DVDT <X DR/VOL(K1) 


(with the proper attention paid to indices) in both Q) and Q,: 


2. To handle PTOT and Qi» we write (3) as 


a 17 p ~ WwW ] 
Q, yi-67 Pe WB © DDDV 


l_ ry ; Oh \ 7 
- ¥1. 67(P + €,Q) ' €4Q5) P WB» DDDV 


67 


(9) 


10) 


(11) 


og 


2 
Q\ 


or 


1 2 : 
asl - 1. 67p 19) - 1. 67(P + £59.) p = (0 (14) 


(WB * DDDV) 


so that 


9 4 1. 67p(P + €5Q,) 
1.67) €, + (1-67) €,)) + ——__5,* 
(WB * DDDV) 


2/(WB » DDDV)2 
let Ee Se ee ge ae 


Q, : (15) 


where the negative-sign solution is discarded because Q, must be 


positive. 


Thus, we should first compute Q, from (4) in which DDDV is 
given by (10), then compute Q, from (15), and finally Q from (2). 


oe 569 Initial Total Energy After Rezone 


If the rezone flag (IENCHK) equals 1, Subroutine REZONE calls 
Subroutine ENECHK to compute a new value for the total energy in the 


rezoned column. 


3.5.10 Resetting Vertical-Direction Indices of Time-Dependent 


Event Points 


After the basic portion of the rezoning has been completed, the 
index (KX(L)) for the cell, if any, containing the event-point L is com- 


puted, as well as the fractional-distance location (FRC(L)) within the 
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cell. These quantities, KX(L) and FRC(L), are passed through EVENTX 
Common. If the point of interest is above the rezoned column, a message 


to that effect is printed. 


3.5.11 Hydro Cycle Counters 


At the end of Subroutine REZONE, provided at least one column 
has been rezoned, the status of the hydro cycle counters (KCYC(NN)) is 


printed and then the counters are zeroed. 


3.6 VERTICAL MOTION OF A POINT (SUBROUTINE TIMVAR) 


To meet the needs of PDI, in their overlay fireball model, to 
follow the motion of a specified Lagrangian point (in particular, an 
event-related point), we developed Subroutine TIMVAR. This routine 
computes the time-dependent position of a Point X and its corresponding 
velocity at the evaluation time. In application, we have followed only 
the event points, but the program could be easily modified to follow the 
motion of any desired points, e.g., the positions of selected debris 


locations. 


As input, Subroutine TIMVAR needs the three indices (IX, JX, 
KX) of the cell containing the point to be followed and the fractional 
vertical-position (FRC) of the point within the cell. The quantities KX 
and FRC are reset at rezone. The position and velocity of the point at 


an evaluation time are linearly interpolated from the positions and 


velocities of the boundaries of the cell containing the specified point. 


4. NRLHYD (DIFFERENTIAL-QUADRATURE FORM OF 
LAGRANGIAN HYDRODYNAMICS) 


4.1 INITIALIZATION OF LAGRANGIAN POINTS IN A COLUMN 
UNDER AMBIENT CONDITIONS 
The initialization of the Lagrangian points in a column is started 
by a call from Subroutine GRIDON to (the NRLHYD form of) Subroutine 
ATMOSG. 


4.1.1 Lagrangian- Point Altitude 


We have placed a NRLHYD Lagrangian point at the center of 
each of the SAIHYD cells plus a (permanently) stationary point at the bot- 
tom boundary (HABOT) of the SAIHYD column of cells. Thus, we must 
reset the value of NCELZ used in SAIHYD according to the expression 


(NCELZ) ver Hyp = (NCELZ) 64 typ +1 


To provide for 'cell boundaries’ that can be used for energy-deposition 
purposes in NRLHYD, we reset the value of HABOT used in SAIHYD ac- 
cording to the expression 


0.25 * DELHA(1) 


(HA BOT) = (HABOT) 


NRLHYD SAIHYD ~ 


4.1.2 Lagrangian- Point Velocity 


The velocity of each of the Lagrangian points is set to zero. 


4.1.3 Lagrangian- Point Pressure 

The pressure at each of the Lagrangian points is set by the 
value of the ambient pressure (PAMB) passed through ATMOUP Common 
by a call to Subroutine ATMOSU(2, HCKM), with HCKM the altitude of 


the Lagrangian point. 


4.1.4 Lagrangian- Point Density 


The mass density at each of the Lagrangian points is set by the 
value of the ambient density (RHOAMB) passed through ATMOUP Common 
by the call to Subroutine ATMOSU(2, HCKM). 


4425 Lagrangian- Point Species Densities and Other Chemistry 
Quantities 


Species densities and other chemistry quantities are not required 
for the hydrodynamics but they are for the chemistry. In general, these 
quantities are determined in Subroutine ATMOSG at the beginning of a 
loop over the cells in a column by three calls to Subroutine ATMOSU(2, 
HCKM), I[ONOSU(2, HCKM), and SPCMIN(2, HCKM), where HCKM is the 
altitude (in kilometers) of the Lagrangian point. We list in Table 1 (in 


Section 3.1.7) the chemistry quantities and how they are determined. 


4.1.6 Weighting Factors for Local Truncation Error in NRLHYD 
Integrator 
Weighting factors, called HEVEPS(K) in Subroutine ATMOSG, 
are set for estimating local truncation errors in the NRLHYD integrator 
(Subroutine RKAM2G), where they are called AUX(K, 6). For integrat- 
ing the acceleration to obtain the velocity, the weighting factors are 


100 1 


HEVEPS(K) T01l NCELZ, 


For integrating the velocity to obtain the radius, the weighting factors 


are 


HEVEPS(K) rai NOEL 


4.2 CONTROL ROUTINE FOR ADVANCING THE 
HYDRODYNAMICS (SUBROUTINE HYDROG) 


Subroutine HYDROG, the routine controlling the subroutines 


doing the actual advancement of the hydrodynamics, performs the fol- 


lowing functions: 


a. 


Establishes the proper headings for the two sets of time- 
slot data by testing on the value of the parameter (IDREG) 
for the number of hydro and chemistry updates. IDREG, 
initially set and updated in DRIVER, is passed through 
PRPREG Common. 


Controls the looping over all the columns in the grid while 
the hydrodynamics is being done. 


Zeroes the emergency-rezone flag (IEMREZ). 


Effects the transfer of the NCOL Lagrangian-point quanti- 
ties for the column of interest from large core memory 
(LCM) to small core memory (SCM). This transfer is 
done by a call to Subroutine ECRD(BUF, IWI, NCOL), 
where 


IWI = (NN - 1) ¥ NCOL+NTIME1 , 


so that NCOL quantities starting at index IWI in LCM are 
transferred to the BUF array in SCM. 


Initializes the following variables before calling Subroutine 
INLHVE which initializes the variables for the NRLHYD 
integrator (Subroutine RKAM2G): 


Geocentric radius of Lagrangian-point K: 


Y(N + K) = BUF(LOC + 1) + RE 
ZZ(K) Y(N + K) 


iH] 


Velocity of Lagrangian-point K: 


V(K) = BUF(LOC + 5) 
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Weighting factors used for local truncation error in Sub- 
routine RKAM2G: 


HEVTEM(K) = HEVEPS(K) 
HEVTEM(N + K) = HEVEPS(N + K) 


Mass density of Lagrangian-point K: 
RHOO(K) = BUF(LOC + 2) 

Pressure at Lagrangian-point K: 
PO(K) = BUF(LOC + 15) 


Calls Subroutine INLHVE to initialize the differential quad- 
rature variables for the column of interest. 


Calls Subroutine RKAM2G to advance the hydrodynamics 
in a column from the start time (TIME) to the stop time 
(TSTOP). 


Provides a procedure for an emergency rezone of the 
column, in an attempt to keep the problem running after 11 
halvings of the time step in Subroutine RKAM2G, indicated 
by IHLF = 11, fail to enable the integrator to achieve the 
desired accuracy. 


Prints a message if one emergency rezone is insufficient 
for Subroutine RKAM2CG to achieve the desired accuracy. 


Checks whether or not Column NN needs to be rezoned, 
based on the value of the geocentric radius of the topmost 
Lagrangian point relative to the criterion radius (HGTMAX), 
and sets the rezone flag (IREZ(NN)) if required. ITREZ(NN) 
is incremented by one if rezoning is required for Column 
NN, thus making its value odd. After rezone is performed 
in Subroutine REZONE, based on IREZ(NN) being odd, 
IREZ(NN) is then incremented by one in Subroutine REZONE, 
thus making its value even. The number of rezones already 
performed on Column NN can be found by dividing IREZ 

by two. 


EE 


h.2. Prints a message, if a rezone is required, and states the 
rezone number and the associated time. 


i. Outputs the following hydrodynamics quantities for each of 
the Lagrangian points in Column NN at time TSTOP: 


(1) Index 

(2) Altitude 

(3) Mass density 

(4} velocity 

(5) Retio of current density to initial density 
(6) Specific internal energy 

(7) Pressure 


j;. Updates the BUF array (compare Step e. 1 above). 
k. Effects the transfer of the NCOL Lagrangian-point quanti- 
ties for the column of interest from SCM to LCM by a call 


to Entry Point ECWR(BUF, IWI, NCOL) in Subroutine ECRD. 


1. Sets and prints the values of the altitudes and velocities of 
the time-dependent event points by calling Subroutine 


TIMVAR. 
4.3 INITIALIZATION ROUTINE FOR NRLHYD (SUBROUTINE 
INLHVE) 


Subroutine INLHVE initializes the variables in the differential- 
quadrature treatment of the Lagrangian hydrodynamics by computing at 


each of the Lacrangian points in the column the initial values of 
a. The density ratio, RRO(K) = 1. 
b. The pressure-to-density ratio, CSQO(K) = PO(K)/RHOO(K). 
c. The (negative) logarithmic pressure derivative, 


N 
DPDZO(K) = - )> COEF(K, L) = In[PO(L)] , 
L=1 


where COEF(K,L) is the array of coefficients in the 
Lagrange differentiation formula, set by a call from 


te. 


—- a 


Subroutine INLHVE to Subroutine SETFUN and passed 
through VAR1 Common. 


4.4 ROUTINE FOR LAGRANGE-DIF FERENTIATION- FORMULA 

COEFFICIENTS (SUBROUTINE SET FUN) 

Subroutine SETFUN(M) computes the coefficients COEF(I, J) in 
the Lagrange differentiation formula, where M is the maximum num- 
ber of points used in calculating the coefficients. A value of five for 
M is used for the third and higher-number of points from either end of 
the column of points. For the second point from an end, a three-point 
formula is used. For an end point, a one-sided two-point formula is 
used. The coefficients in Row I of the COEF(I, J) array are used to 
calculate the derivative of a function at Lagrangian-point I in the 


column. 


The Lagrange interpolation formula for a polynomial of degree 


m-1 [AS-64, p. 878| is 


m 


f(x) dic, ) {(x,) 
j=1 . 


where 


m, 
(Xx - X)) 
k=1 k 
mn, a 
(x, - x) 


eet eee ree a gn ee ee . eee 


The Lagrange differentiation formula |AS-64, p. 882], obtained 


by differentiating the Lagrange interpolation formula, is 


m 
df , 7 
— f(x) = : 1 Oe) 
ce ) de f(s) (x;) 
j=1 
where 
m_, 
(x - x,) 
m , = 
c.(x) = Pe ma Mm, 
L=1 
ij (x - x,) : (x. - x,) 
kAj 


m_, 
a ~ 2) 
i k 
aa. ee 
Z J 
ij C(x) ta] oe 
43 (x. - x,) [| (x, - X;) 
LZi ’ 
kA} 


This formula is given by Gardner |Ga-74b, GP-74, GP~75] in the fol- 


lowing form: 


m_, 
(x; - X,) 
— 4 diet 
s a 5 
“ij . a X,- X, = m_, it 
Lei] (x; - x) 
k=1 
k# jet 


The coefficient at is called COEF(U, J) in Subroutine SETFUN. 
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4.5 INTEGRATOR FOR THE NRLHYD DIFFERENTIAL EQUATIONS 

(SUBROUTINE RKAM2G) 

Subroutine RKAM2G evaluates the velocity and location of the 
Lagrangian points ina column. The integration in time is performed 
by using a Hamming-type modification [LS-71c, p. 155] to the Adams- 
Bashforth two-step second-order predictor equation |LS-71c, p. 180, 
Table 4.1 (A-B, k=2)| and the Adams-Moulton one-step second-order 
corrector equation |LS-71c, p. 180, Table 4.1 (A-M, k=1)|, which 
results in an overall modified predictor-corrector procedure of third 
order |LS-71c, p. 157]. The required two-step starting values for the 
modified predictor-corrector scheme are obtained from a second-order 
Runge-Kutta method |LS-71c, p. 41, Eq. (2.1-11)], supplemented by a 


corrector calculation. 


Subroutine RKAM2G, originally prepared by NRL, has been 


modified by SAI in the following principal respects: 


a. Reduction of the value for the time step (PRMT(3)), set as 
a Data Statement in DRIVER, from 10 sec used in NRLHV4 
to 1 sec used in earlier versions, as found necessary by 
SAI in several test problems to avoid having the Runge- 
Kutta portion of RKAM2G abort due to negative numbers 
being raised to powers in Subroutine FUNCT. 


b. Removal of the constraint on the modified-predictor- 
corrector portion of RKAM2G that it could not double the 
time step with respect to PRMT(3). 


c.1. Removal of the tacit constraint, for proper working, that 
the time interval from the start-time to the stop-time was 
exactly divisible by the time step. 


c.2. Addition of the capability to interpolate to obtain the values 
of the velocity (Y(K)) and position (Y(N+K)) at the stop-time 
(X=PRMT(2)) when the integrator oversteps the specified 
stopping time. 
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Reduction of the value for the bound on the truncation error 
(EPSSET = PRMT(4)), set as a Data Statement in DRIVER, 
from 1000 cm used in NRLHV4 to 100 cm used in earlier 
versions, to help prevent aborting. 


e. Addition of a hydro-cycle counter. 


f. Elimination of the capability to integrate backwards in 
time. 


g. Addition of a test in Subroutine HYDROG to detect when 
Subroutine RKAM2G, in either the RK or AM2G portion of 
the routine, fails to achieve the desired integration ac- 
curacy after 11 halvings of the time step and returns 
clandestinely to Subroutine HYDROG. In sucha circum- 
stance, after the proper saving and resetting of rezone 
flags, Subroutine HYDROG calls Subroutine RE ZONE to 
perform an emergency rezone on the column to help relieve 
the awkward conditions facing the integrator. In this re- 
zone the top and bottom Lagrangian points are kept fixed 
and the interior points are uniformly spaced. To achieve 
the emergency rezone, a proper shuffling of the column 
data is made to Time-Slot 2 for starting the rezone and 
back to Time-Slot 1 after the rezone for later use by the 
chemistry routines. Subroutine RKAM2G is then permitted 
to integrate for the column. Only one such rezone at a 
eiven time is potentially useful, so an exit is called if the 
program attempts to repeat such an emergency rezone. 


h. Revision, and extensive addition, of comment cards. 


4.6 DIFFERENTIAL EQUATIONS FOR HYDRODYNAMICS IN 
NRLHYD (THE DERIVATIVE ROUTINE, SUBROUTINE FUNCT) 
Subroutine FUNCT computes the derivatives (the velocity and 

acceleration at the Lagrangian points) for the ordinary differential equa- 

} tions (for the position and velocity at the Lagrangian points) integrated 

in Subroutine RKAM2G. The ordinary differential equations have been 
obtained |Ga-74a, Ga-74b, GP-74, GP-75| from the partial differential 


equations of hydrodynamics by assuming isentropic (shock-free) flow and 


by replacing the partial derivatives in the spatial direction by derivatives 
of en analytic fit function, in accordance with the method of differential 
quadrature |BK-72b|. Asa result of these assumptions, the acceleration 


at a Lagrangian point can be expressed in terms of 


a. The initial values of 

(1) The pressure-to-density ratio at the point and 

(2) The logarithmic pressure derivative at the point and 
b. The current values of 


(1) The ratio of the current Lagrangian geocentric radius 
to the initial Lagrangian geocentric radius, 

(2) The ratio of the current density at the point to the 
initial density at the point, and 

(3) The first and second (ordinary) derivatives of the cur- 
rent Lagrangian geocentric radii with respect to the 
initial Lagrangian geocentric radii, evaluated by dif- 
ferentiating polynomial fit functions for the current 
Lagrangian geocentric radii. 


The equations of motion in NRLHYD, as written by Gardner 
|\Ga-74b, GP-74, GP-75]| ina form requiring only the initial pressure, 


density, and velocity at the initial fluid element locations, are 


FAS, v(r, t) 
at 
av(r,t) _ _ Pol”? [Ri rt)? perce, ty 
ot Prt r or 
open aid (eRe! o*Rir,t)  2y | Rr Uy oR(r,t) _ | 
or a Gr by" a ed or { 


1 initial geocentric radius of Lagrangian-point K. 
R(r, t) current geocentric radius of Lagrangian-point K. 


v(r,t) current velocity of Lagrangian-point K. 
current acceleration of Lagrangian-point K. 
p (r) = initial pressure at Lagrangian-point K. 

o (i) initial density at Lagrangian-point K. 


The correspondence between these variables and the Fortran variables 
is given in Table 3, where we have also given the expressions used to 

evaluate the variables; of particular importance is the evaluation of the 
spatial derivatives by use of the polynomial fit functions, in accordance 


with the differential-quadrature method. 
In terms of the Fortran variables, the equations of motion 
would be expressed as follows: 
DERY(N + K) Y(K) 
aa | 


DERY(K) CSQ0(K) » REO) 5 


(A3) 


DPDZO0(K) 


. GAMMA » [a2 x AY ‘oa ls : i} - GRAV 


However, in the NRL Subroutine FUNCT, there are two depar- 
tures from the description in GP-74 and GP-75. First, the first and 
second derivatives of the Lagrangian positions (for interior points) are 
computed from a weighted-average procedure. In addition to the deriva- 


tives computed by use of the polynomial fit-functions mentioned above, 
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derivatives are also computed by taking first and second differences. 
The contribution of the direct-difference derivative to the derivative 
ultimately used depends on a weighting factor (CHI) so chosen that it 1s 
unity at the starting time and for any compression and decreases from 


unity as expansion occurs. 


Thus, NRL has introduced the following quantities: 


A2 = 1/Max|RR(K), 0. 25] 


. 2 2 
; |... 2eltks 1) - 228) ( 27,(K) - ZZ(K - 1) 
CHy)= Mims wax larg Ke 1)- YING RK)’ VYING K)- YN TK - a) | 


ag -1[_22(K+1)- 22(K)_, _2Z(K) - ZZ(K - 1) 
a(Y(N+K+1)- YIN+ K) Y(N+ K)- Y(IN+ K- 1) 


A2 (1 - CHI) x A2 + CHI x A4 


Y(N + K+ 1) - Y(N + K) 
ZRe(K + 2) = ZAK 


_ ¥(N + K) - Y(N+ K- 1) 2 
ZZ(K)- ZZ(K-1) | ZZ(K+1)- ZZ(K- 1) 


Al (1 - CHI) x Al + CHI X AS 


It is these CHI-weighted expressions for Al and A2 that are actually used 


in the acceleration equation. 


The second departure is that, to partially compensate for the 
simplifying assumption of isentropic flow, NRL has int roduced an arti- 
ficial viscosity term at all interior points to help maintain a smooth 
solution, i.e., on the right-hand side of the acceleration equation there 


is added the term 
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)| DVDZP| x DVDZP - |DVDZM| » DVDZM, ; 


ANUC »* 


where 
ANUC ANU > [ZZ(K t 1) - ZZ(K - 1) | ; RRO) 
ANU Oot 


DERY(N + K+ 1) - DERY(N + K) 


ITIW7P = 
DVDZI Zot = Sk 


and DVDZM equals the value of DVDZP from the previous pass through 


the DO-loop over the points in the column. 


In addition to the above-mentioned departures, SAI has made 
two changes in the NRL Subroutine FUNCT. First, we have modified 
the NRL procedure at the (stationary) bottom Lagrangian point by requir- 
ing the density to remain constant instead of letting it be determined by 
the amount of motion undergone by the second-lowest point. The treat- 
mer of the bottom-point density obviously affects the HAG-interpolated 
values for the density between the two lowest Lagrangian points but (sur- 
prisingly) it does not affect the motion of any of the moving Lagrangian 
points until the column has been either restarted (with a new initializa- 
tion) or rezoned (also with a new initialization). With a new initializa- 
tion there is a new pressure profile, partially determined by the 
boundary condition prescribed for the density at the lowest point. It is 
this new pressure profile that affects the motion of the points above the 


lowest. 


Our second change in Subroutine FUNCT is that we have pro- 
vided a test to determine whether or not RRO(K), the ratio of the current 
density to the initial density at Lagrangian-point K, is negative as a 
result of two adjacent points crossing owing to too large a time step. 


(If RRO(K) is negative, the program used to abort without having a chance 


ee —— - — “ rage ee -- a a - --- — age 
e $$ $$ 


to try a smaller time step.) To prevent an abortion when RRO(K) is 
negative, we require the program to return to near the beginning of 
Subroutine RKAM2G to try again with one-tenth the previous time step. 
Two such emergency returns are permitted before allowing the column 
to undergo an emergency rezone, initiated by artificially setting the 


bisection counter IHLF to 11 and causing a return to Subroutine HYDROG. 


4.7 REZONE IN NRLHYD (SUBROUTINE REZONE) 


ra te | Purpose 


To provide for maintaining a reasonable resolution of the 
Lagrangian points in the (more important) lower portion of the HAG 
module as a column expands to very high altitude, NRL developed a re- 
zone procedure which we modified and incorporated into Subroutine 
REZONE. 


r= 
~] 
a] 


General Features of Subroutine REZONE 


Subroutine REZONE resets the positions and properties of the 
Lagrangian points for each column whose geocentric radius of the upper- 
most point exceeds a given geocentric radius HGTMAX at the hydro- 
advance time TSTOP. The number of Lagrangian points in a rezoned 
column is kept the same (NCELZ yp NCELZ., | 
Lagrangian-point separation distances (DELHA2(K)) are read by DRIVER 
in Namelist SETUP1 and passed through REZONN Common. If 


1). The new 


DELHA2(1) = 0., the rezoned separation distances are evenly Spaced as 
proposed by NRL. Table 4 gives the currently used values for the initial 


and rezoned Lagrangian-point separation distances and altitudes. 


Subroutine REZONE uses linear interpolation for velocity, 


electron temperature, and N, vibration temperature and uses logarithn ¢ 
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interpolation for pressure and density-dependent variables. These in- 


terpolations are made by calling Subroutine INTERP. 


Subroutine REZONE is called from DRIVER after all the col- 
umns have been advanced first for hydrodynamics and then for chemistry. 
In looping over the columns, Subroutine REZONE rezones only those 
columns for which a rezone flag (IREZ(NN)) has been set (based on the 
rezone criterion) in Subroutine HYDROG immediately after Subroutine 


RKAM2G completes the hydrodynamics for Column NN. 


If the rezone flag has been set for Column NN, then Subroutine 
ECRD is called to obtain the column quantities from Time-Slot 2. For 
development purposes, we print the status of the Lagrangian-point quan- 
tities, and also the molecular-ion density and heavy-particle tempera- 
ture, before rezone. After the new Lagrangian-point quantities are 
determined, Entry-Point ECWR in Subroutine ECRD is called to store 
the Lagrangian-point quantities in Time-Slot 2. Again, for development 
purposes, we print the rezoned Lagrangian-point properties, including 


the molecular-ion density and heavy-particle temperature. 


4.7.3 Resetting Vertical-Direction Indices of Time-Dependent 

Event Points 

After the basic portion of the rezoning has been completed, the 
index (KX(L)) for the Lagrangian points, if any, associated with the 
event-point L is computed, as well as the fractional vertical-position 
parameter (FRC(L)). These quantities, KX(L) and FRC(L), are passed 
through EVENTX Common. If the point of interest is above the rezoned 


column, a message to that effect is printed. 
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4.7.4 Hydro Cycle Counters 


At the end of Subroutine REZONE, provided at least one col- 
umn has been rezoned, the status of the hydro cycle counters (KCYC(NN)) 


is printed and then the counters are zeroed. 


4.8 INTERPOLATION ROUTINE IN NRLHYD (SUBROUTINE 

INTERP) 

Subroutine INTERP(K, DUMX, DUMY, M, X, Y, MODE) is 
called from Subroutine REZONE to provide interpolated values of the 
Lagrangian-point quantities. The definitions of the arguments are as 
follows: 

K = Index for which the abscissa value (X), at which one wishes 


io determine the ordinate value (Y), satisfies the relation 
X(K) < X < X(K+1). 


| 
DUMX = Array of abscissas (old Lagrangian-point altitudes). 
DUMY = Array of ordinates (variable dependent on old Lagrangian- 
point altitudes). 
M Dimension of DUMX and DUMY. 
X X- coordinate at which interpolated value for Y is desired 
: (new Lagrangian-point altitudes). 
i Y Returned interpolated value. 
| MODE 1, linear two-point interpolation. 
. 2, logarithmic two-point interpolation. 
i 
f 4.9 VERTICAL MOTION OF A POINT (SUBROUTINE TIMVAR) 
st 
4 See Section 3.6 in SAIHYD for a general discussion of the prop- 


erties of Subroutine TIMVAR, since the NRLHYD version is very similar 
to the SAIHYD version, except for two related features. Instead of 


vertically locating the Point X as having a fractional-distance location 
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within a cell as in SAIHYD, we use the fractional-distance between the 


two closest Lagrangian points. We also specify the index of the closest 
Lagrangian-point K, with an assigned plus or minus sign to denote 


whether the Point X is above or below Lagrangian-point K, respectively. 
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5. CONTROL ROUTINE FOR ADVANCING THE CHEMISTRY | 
(SUBROUTINE CHEMG) | 


Subroutine CHEMG, the routine controlling the subroutines 


doing the actual advancement of the chemistry, performs the following 


functions. 
a. 
ae 
Cc. 
d. 
RHO(1) 
TEL(1) 
ENE(1) 
EN4S(1) 
EN2D(1) 
OATOM(1) - 
EN2(1) 
O2(1) 
ENO(1) 
ENP(1) 
OP(1) 


PRES(1) 
TVIBN2(1) 
CO2(1) 


Controls the looping over all the columns in the grid while 
the chemistry is being done. 


Effects the transfer of the NCOL cell quantities for the 
column of interest from large core memory to small core 
memory, Similar to the transfer effected from HYDROG in 
preparation for doing the hydrodynamics. In addition to 
transferring the just-described data from Time-Slot 1, 
data are also transferred from Time-Slot 2 since some of 
these data are required for the GET-chemistry module. 


Prints, for development purposes, the status of the 
Lagrangian-cell quantities, and also the molecular-ion 
density and heavy-particle temperature, before doing the 
chemistry. 


Initializes the variables in SPECEF Common before calling 
Subroutine CHEMEF. These variables for Cell K 
(Kk = 2, NCELZ) are: 


BUF1(LOC+ 2) Density at Time T1 
BUF1(LOC+ 3) Electron temperature at Time T1 
BUF1(LOC+ 4) Electron density at Time T1 
BUF1(LOC+ 7) N(4s) density at Time T1 
BUF1(LOC+ 8) N(2D) density at Time T1 
BUF1(LOC+ 9) O density at Time T1 
BUFI1(LOC+10) Nog density at Time T1 

- BUFI(LOC+11) Og density at Time T1 
BUF1(LOC+12) NO density at Time T1 
BUF1(LOC+13) N* density at Time T1 
BUF1(LOC+14) O* density at Time T1 


BUF1(LOC+15) 
BUF1(LOC 
BUF1(LOC+18) 


Pressure at Time T1 
N2 vibration temperature at Time T1 
CQ 9 density at Time T1 
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QAM B(1) 
HE (1) 
RHO(2) 
PRES(2) 


We note: 


BUF1(LOC+20) Ambient ion production rate at Time T1 
BUF1(LOC+22) He density at Time T1 

BUF2(LOC+ 2) Density at Time T2 

BUF2(LOC+15) Pressure at Time T2 


(1) He is not used in the GET-chemistry but is used in the 


NRL-chemistry module, (2) the electron temperature and N@2 vibrational 
temperature are the same quantity in the GET-chemistry but are differ- 
ent in the NRL-chemistry module, and (3) pressure at Time T2 is used 

in the GET-chemistry but is not used in the NRL-chemistry module. 


e. 


Calls Subroutine CHEMEF to advance the chemistry ina 
column from the start time T1 to the stop time T2. 


Computes, for the benefit of Subroutine ENECHK in SAIHYD 
if the energy-check flag (IENCHK) equals one, the total 
energy returned to the cell by the chemistry. 


Updates the BUF2 array (compare Step d above). 


Outputs the following chemistry quantities for Column NN 

at Time T2: Electron temperature; N92 vibrational tempera- 
ture; densities of electrons, N(4s), N(2D), O, Na, O92, NO, 
N*, and O*; pressure; heavy-particle temperatures at 
Times Tl and T2. 


Effects the transfer of the NCOL cell quantities for the col- 


umn of interest from small core memory to large core 
memory by a call to Entry Point ECWR in Subroutine ECRD. 
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6. TEST PROBLEM AND RUNNING TIMES 


To compare the two hydrodynamics modules, SAIHYD and 
NRLHYD, we present results for a three-minute problem in which there 
is a large-energy event at 200-km altitude at t = 0 sec. The grid was 
limited to a row of six columns, spaced 100-km apart at 100-km altitude. 
The atmosphere used corresponds to the CIRA-1965 atmosphere for 
which the solar flux, SBAR, equals 157.57 (see Vol. 14a). The 
event point was placed at the center of Column 6 (achieved by resetting 
the event point after the grid was established). The problem was run 
with both SAIHYD and NRLHYD, but with only one chemistry module, 
GETCHM. 


Breakdowns of the running time and the hydrodynamic cycles 
for this problem are given in Tables 5 and 6, respectively. SAIHYD was 
about 32 times faster than NRLHYD for this problem. For the two-event 
sample problem presented in Vol. 13, SAIHYD was about four times 


faster than NRLHYD. 


Figures 13 and 14 show the altitude-vs-time profiles of every 
other Lagrangian cell-center (or point) for Columns 6 and 3, respectively, 
for both SAIHYD and NRLHYD. The motions predicted by the two forms 


of the hydrodynamics seem to agree quite well, on the whole. 


Figure 15 shows the electron density profile in Column 6 at 
t = 60 sec for both SAIHYD and NRLHYD, and Fig. 16 shows the same 
information for Column 38 at t = 120 sec. The constancy of the electron 
density over a given altitude interval corresponds to the assumed uni 


formity of the cell properties throughout a cell. 


Table 5. Breakdown of SAI-HAG-Module Running Time on CDC-7600 
for One-Event Test Problem for SAIHYD and NRLHYD.® 
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At, sec 


0-1 
1-3 
3-10 
10-30 
30-60 
60-90 
§0-120 
120-150 
150-180 


0-1 
1-3 
3-10 
10-30 
30-60 
60-90 
90-120 
120-150 
150-180 
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Table 6. 


Cycles for the Two Versions of the Hydrodynamics Module 


Used in the One-Event Test Problem. 


A Cycles for Column 
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Cumulative Cycles for Column 


SATHYD 
6 1 
2 2 
3 3 
8 5 
19 9 
18 12 
5 15 
4 18 
B) 21 
4 24 
NRLHYD 
4 2 
8 6 
26 14 
63 25 
o4 39 
19 53 
19 67 
14 81 
14 95 
93 
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Column 6 
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13. Altitude-versus-Time Profiles of Lagrangian Positions in 
Column 6 of Test Problem for SAIHYD and NRLHYD. 
K labels the vertical cell number in SAIHYD; only odd 
values are shown. 
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Column 3 


1100 o SAIHYD 
o NRLHYD 
1000 e SAIHYD & NRLHYD 


Altitude, km 


SAIHYD Rezone NRLHYD Rezone 


0 30 60 90 120 150 180 
Time, sec 


Fig. 14. Altitude-versus-Time Profiles of Lagrangian Positions in 
Column 3 of Test Problem for SAIHYD and NRLHYD. 
K labels the vertical cell number in SAIHYD; oniy odd 
values are shown. 
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= 120 sec 


Column 3 


—— SAIHYD (before rezone) 
-—-—NRLHYD 


‘---r ere 


900 


£700 
o 
'S 
= 
a 500 

300 

100 

104 10° 10° 10! 10° 
Electron Density, em73 
Fig. 16. Electron Density Profiles in Column 3 at t = 120 sec of 


Test Problem for SAJHYD and NRLHYD. 


For the reader who desires more detailed results for the hy- 
drodynamics (and chemistry), we are providing the computer printouts 
for both the SAIHYD and NRLHYD runs. Before doing so, we provide 
the following general guide. At t = 0. 0 , before the energy deposition 
from the event, we print the cell quantities for ambient (gridded) at- 
mospheric conditions. A guide to this print, as to the input to the hydro 
and chemistry at t = 0. 0° and to the output from rezone, is given in 


Table 7. The entries in the blocks of hydro and chemistry quantities 


are defined in Tables 8 and 9. The computer printout is given in 


Tables 10 and 11 for SAIHYD and NRLHYD, respectively. 


Table 7. Guide to Printout for Cell Quantities. 


No. Heading _ 


eee ae Quantity _ 
1 HC Cell-center altitude, cm |[SAIHYD | 
HC Altitude of point, cm |NRLHYD] 
2 RHO Mass density, g/cm? 
3 T ELECT Electron temperature,” K 
4 ENE Electron density, on 
5 VCENTR Cell-center velocity, cm/sec [SAITHYD | 
VPOINT Velocity of point, cm/sec [NRLHYD] 
6 VEL Cell-top-boundary velocity, cm/sec [SAIHYD] 
-- Not used [NRLHYD] 
7  EN4S Density N(4s), em7? 
8 EN2D Density N(“D), em7° 
9 O Density O, em7? 
10 N2 Density No» cm” 
11 02 Density Oo, cm 3 
12 ENO Density NO, cm 
13. ENP Density N*, em7? 
14 OP Density Or, oe" 
is P Total pressure, dyne/em? 
16 Q Artificial viscous pressure, dyne/cm2 [SATHYD] 
ne Not used [NRLHYD] 
17 T N2VIB = T ELECT 
18 CO2 Density CO,, em** 
19 -- Not used 
20 QAMB Effective ambient ion-pair production rate, (ion-pair)/em™? 


21 (HTOP)O Altitude of cell-top-boundary in ambient atmosphere, cm [SAIHYD] 
-- Not used [NRLHYD] 

22 HE Density of He, en? 

23 M’ Density of molecular ions (e-N*-O*), om”? 


24 TH Heavy-particle kinetic temperature, “K 


E T 3 ; ; 
* Actually, temperature common to electrons, O( D)-to-O(’ P) population ratio, and Ny 
vibration. 
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Table 8. 
Column 
No. Heading 
1 K 
2 CEN. ALT. 
ALTITUDE 
3 DENSITY 
4 VCENTR 
VEL 
5 VOLUME 
DEN. RATIO 
6 INT. E. 
7 PRESS. 
8 Q 
9 TOP ALT. 
10 VTOP 


Guide to Printout for Hydrodynamic Output Quantities. 


Quantity 
Cell index 
Cell-center altitude, cm [SAIHYD] 
Altitude of point, cm [NRLHYD] 
Density, e/em* 
Cell-center vertical velocity,” cm/sec [SAIHYD] 
Velocity of point, cm/sec [NRLHYD] 
Steradianal volume of cell, em/sr [SAIHYD] 
Ratio of current- to initial-density [NRLHYD ] 
Internal energy, erg/g 
Pressure, dyne/cm 
Artificial viscous pressure, dyne/cm2 [SATHYD | 
Not used [NRLHYD ] 
Cell-top-boundary altitude, cm [SAIHYD] 
Not used [NRLHYD] 
Cell-top-boundary velocity, cm/sec [SAIHYD ] 


Not used [NRLHYD] 


@Mean of two adjacent boundary velocities. Positive sign is upward motion. 
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Table 9. Guide to Printout for Chemistry Output Quantities’ 


(from CHEMEF®). 


Column 
No. Heading _ Quantity 
1 K Cell index 
2 TEL Electron temperature,” °K 
3 TVIBN2 = TEL 
4 ENE Electron density, oats 
5 N(4S) Density N(4S), em7? 
6 N(2D) Density N(2D), em™ 
7 oO Density O, en”? 
8 N2 Density No, on” 
9 02 Density Op, en 
10 NO Density NO, wn 
11 NP Density N’, an 
12 OP Density oO", oa 
13 PRES Pressure, dyne/em” 


14 TEM(1) Heavy-particle temperature at TDAT1,° °K 


15 TEM(2) Heavy-particle temperature at TDAT2,° °K 


4411 quantities are at time TDAT2 except TEM(1) which is at TDAT1. 


bemperature common to electrons, o(p)-to-oPP) population 
ratio, and No vibration. 


“Temperature inferred by subtracting electron pressure from total 
pressure. 


dhe molecular ion density, regarded at [NO*], may be computed 
from charge conservation: [NO*] = ENE - NP - OP. 


“Chemistry developed by GET personnel. 
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